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Abstract 

Using Planck maps of six regions of low Galactic dust emission with a total area of about 140 deg^, we determine the angular power spectra of 
cosmic infrared background (CIB) anisotropies from multipole I = 200 to ^ = 2000 at 217, 353, 545 and 857 GHz. We use 21 -cm observations of 
Hi as a tracer of thermal dust emission to reduce the already low level of Galactic dust emission and use the 143 GHz Planck maps in these fields 
to clean out cosmic microwave background anisotropies. Both of these cleaning processes are necessary to avoid significant contamination of the 
CIB signal. We measure correlated CIB structure across frequencies. As expected, the correlation decreases with increasing frequency separation, 
because the contribution of high-redshift galaxies to CIB anisotropies increases with wavelengths. We find no significant difi'erence between the 
frequency spectrum of the CIB anisotropies and the CIB mean, with A///=15% from 217 to 857 GHz. In terms of clustering properties, the 
Planck data alone rule out the linear scale- and redshift-independent bias model. Non-linear corrections are significant. Consequently, we develop 
an alternative model that couples a dusty galaxy, parametric evolution model with a simple halo-model approach. It provides an excellent fit to 
the measured anisotropy angular power spectra and suggests that a diff'erent halo occupation distribution is required at each frequency, which 
is consistent with our expectation that each frequency is dominated by contributions from diff'erent redshifts. In our best-fit model, half of the 
anisotropy power at ^=2000 comes from redshifts z < 0.8 at 857 GHz and z < 1.5 at 545 GHz, while about 90% come from redshifts z >2 at 353 
and 217 GHz, respectively. 

Key words. Cosmology: observations 



1. Introduction 

In addition to instrument noise, deep cosmological surveys in 
the far-infrared to millimeter spectral range are limited in depth 
by confusion from extragalactic sources (e.g. Blain et al. 1998; 
Lagache et al. 2003; Dole et al. 2004; Fernandez-Conde et al. 
2008; Nguyen et al. 2010). This limitation arises from the high 
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density of faint, distant galaxies that produce signal fluctuations 
within the telescope beam. As a consequence, the cosmic in- 
frared background (CIB), which records much of the radiant 
energy released by processes of structure formation that have 
occurred since the decoupling of matter and radiation follow- 
ing the Big Bang (Puget et al. 1996; Hauser & Dwek 2001; 
Dole et al. 2006), is barely resolved into its constituents. Indeed, 
less than 10% of the CIB is resolved by Spitzer at 160 pm 
(Bethermin et al. 2010a), about 10% by Herschel at 350 pm 
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(Oliver et al. 2010) and a negligible fraction is resolved by 
Planck^ (Femandez-Conde et al. 2008). Thus, in the absence of 
foreground (Galactic dust) and background (cosmic microwave 
background, CMB) emissions, and when the instrument noise is 
subdominant, maps of the diffuse emission at the angular reso- 
lution probed by the current surveys reveal a web of structures, 
characteristic of CIB anisotropics. With the advent of large area 
far-infrared to millimeter surveys (Herschel, Planck, SPT and 
ACT), CIB anisotropics constitute a new tool for structure for- 
mation and evolution study. 

Cosmic infrared background anisotropics are expected to 
trace large-scale structures and probe the clustering properties of 
galaxies, which in turn are linked to those of their hosting dark 
matter halos. Because the clustering of dark matter is well under- 
stood, observations of anisotropics in the CIB constrain the rela- 
tionship between dusty, star-forming galaxies and the dark mat- 
ter distribution. The connection between a population of galax- 
ies and dark matter halos can be described by its halo occupation 
distribution (HOD; Peacock & Smith 2000; Seljak 2000; Benson 
et al. 2000; White et al. 2001; Berlind & Weinberg 2002; Cooray 
& Sheth 2002), which specifies the probability distribution of the 
number of objects with a given property (e.g., luminosity, stellar 
mass, or star-formation rate) within a dark matter halo of a given 
mass and their radial distribution within the halo. The HOD and 
the halo model provide a powerful theoretical framework for de- 
scribing the connection between galaxies and dark matter halos. 
Once decisions are made about which properties of the halos and 
their environment the HOD depends upon, what the moments of 
the HOD are and what the radial profile of objects within halos 
is, the halo model can be used to predict any clustering-related 
observable. In particular, the halo model predicts that the bias, 
describing the clustering of galaxies in relation to the dark mat- 
ter, becomes scale-independent at large scales. This assumption 
of a scale-independent bias is often made in modelling the CIB. 

The way galaxies populate dark matter halos is not the 
only ingredient that enters into the CIB anisotropy modelling. 
Correlated anisotropics also depend on the mean emissivity 
per comoving unit volume of dusty, star-forming galaxies, that 
results from dusty galaxies evolution models. Such models are 
more and more constrained thanks to the increasing number of 
observations (mainly galaxies number counts and luminosity 
functions), but remain largely empirical. So far, CIB anisotropy 
models have combined (i) a scale-independent bias clustering 
with a very simple emissivity model based on the CIB mean 
(Knox et al. 2001; Hall et al. 2010) or an empirical model of 
dusty galaxy evolution (Lagache et al. 2007) or the predictions 
of the physical model by Granato et al. (2004) for the formation 
and evolution of spheroidal galaxies (Negrello et al. 2007); (ii) 
a HOD with the Lagache et al. (2003) dusty galaxies evolution 
model (Amblard & Cooray 2007; Viero et al. 2009); and (iii) 
a merger model of dark matter halos with a very simple dust 
evolution model (Righi et al. 2008). 

The angular power spectrum of CIB anisotropics has two 
contributions, a white-noise component caused by shot noise 
and an additional component caused by spatial correlations be- 
tween the sources of the CIB. Correlated CIB anisotropics have 



^ Planck (http://www.esa.int/P/flncA:) is a project of the European 
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sortium led and funded by Denmark. 



been measured at 3330 GHz by AKARI (Matsuura et al. 2011), 
3000 GHz by IRAS/IRIS (Penin et al. 2011b), 1875 GHz by 
Spitzer (Lagache et al. 2007; Grossan & Smoot 2007), 1200, 
857, 600 GHz by BLAST and SPIRE (Viero et al. 2009; 
Amblard et al. 201 1), and 220 GHz by SPT (Hall et al. 2010) and 
ACT (Dunkley et al. 201 1). Depending on the frequency, the an- 
gular resolution and size of the survey these measurements can 
probe two different clustering regimes. On small angular scales 
> 2000), they measure the clustering within a single dark mat- 
ter halo and accordingly the physics governing how dusty, star- 
forming galaxies form within a halo. On large angular scales, 
CIB anisotropics measure clustering between galaxies in differ- 
ent dark matter halos. These measurements primarily constrain 
the large-scale, linear bias, b, of dusty galaxies, which is usually 
assumed to be scale-independent over the relevant range. Given 
their limited dynamic range in scale, current measurements are 
equally consistent with an HOD model, a power-law correlation 
function or a scale-independent, linear bias. All models return a 
value for the large-scale bias that is 2-4 times higher than that 
measured for local, dusty, star- forming galaxies (where b ^ 1). 

Owing to its frequency coverage from 100 to 857 GHz, the 
HFI instrument on-board Planck is ideally suited to probe the 
dark matter - star-formation connection. Planck (Tauber et al. 
2010; Planck Collaboration 2011a) is the third-generation space 
mission to measure the anisotropy of the cosmic microwave 
background (CMB). It observes the sky in nine frequency bands 
covering 30-857 GHz with high sensitivity and angular reso- 
lution from 3V to 5'. The Low Frequency Instrument (LFI; 
Mandolesi et al. 2010; BersanelH et al. 2010; Mennella et al. 
2011) covers the 30, 44, and 70 GHz bands with amplifiers 
cooled to 20 K. The High Frequency Instrument (HFI; Lamarre 
et al. 2010; Planck HFI Core Team 2011a) covers the 100, 143, 
217, 353, 545, and 857 GHz bands with bolometers cooled to 
0. 1 K. Polarization is measured in all but the highest two bands 
(Leahy et al. 2010; Rosset et al. 2010). A combination of radia- 
tive cooling and three mechanical coolers produces the temper- 
atures needed for the detectors and optics (Planck Collaboration 
201 lb). Two data processing centres (DPCs) check and calibrate 
the data and make maps of the sky (Planck HFI Core Team 
2011a; Zacchei et al. 2011). Planck's sensitivity, angular reso- 
lution, and frequency coverage make it a powerful instrument 
for Galactic and extragalactic astrophysics as well as cosmol- 
ogy. Early results are given in Planck Collaboration (201 la-u). 

The primary objective of this paper is to measure with 
Planck HFI the CIB anisotropics caused by the clustering of 
star-forming galaxies. To achieve this, we analyse small regions 
of sky with a total area of about 140 deg^, where we are able 
to cleanly separate the foreground (Galactic cirrus) and back- 
ground (CMB) components from the signal. Unlike previous 
CIB anisotropy studies (but see Penin et al. 2011b), we do 
not remove the cirrus by fitting a power-law power spectra at 
large scales, but use an independent, external tracer of diffuse 
dust emission (the Hi gas). We accurately measure the instru- 
mental contributions (noise, beam) to the power spectra of CIB 
anisotropics and use a dedicated optimal method to measure 
power spectra (Ponthieu et al. 2011). All these steps allow us to 
recover for the first time the power spectra of CIB anisotropics 
from 200 < £ < 2000 at four frequencies simultaneously: 217, 
353, 545 and 857 GHz. 

The paper is organized as follows: in Sect. 2 we present 
the data we are using, the field selection and the removal of 
foreground and background components (CMB, Galactic cirrus, 
bright point sources) from the CIB. In Sect. 3 we discuss the dif- 
ferent contributions to the power spectra of the residual maps. 
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Table 1. CIB field description: centre (Galactic coordinates), size, mean and dispersion of Hi column density. 






Figure 1. From left to right and top to bottom: N 1 , AG , SP , LH2 and boot e s fields overlaid on IRIS 1 00 yum map (Miville-Deschenes & Lagache 
2005). Fields Bootes 1 and 2 are both included in the large rectangle. All IRIS images have the same dynamic range, with a linear colour scale 
ranging from dark red to white from to 2 MJy sr~^ . 
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Section 4 describes how we estimated the power spectrum, its 
bias, and errors. Our main results are presented in Sect. 5. This 
section also describes our modelling and discusses the cluster- 
ing of high-redshift, dusty galaxies. We conclude in Sect. 6. 
In the appendices we show two flow charts summarizing the 
data processing and cleaning, and the power spectra measure- 
ments (Appendix A), and we give some details about the dusty 
star-forming galaxy evolution model (Appendix B) and the halo 
model (Appendix C) we are using. Throughout the paper we use 
the WMAP7 cosmological parameters for standard ACDM cos- 
mology (Larson et al. 201 1). 

2. Selected fields and data cleaning 

2.1. Planck data 

We used Planck channel maps of the HFI at 5 frequencies: 
143, 217, 353, 545 and 857 GHz. Their characteristics and 
how they were created is described in detail in the companion 
paper on HFI early processing (hereafter HEP; Planck HFI 
Core Team 2011a). In summary, the channel maps correspond 
to temperature observations for the two first sky surveys by 
Planck. The data are organized as time-ordered information, 
hereafter TOI. The attitude of the satellite as a function of time 
is provided by two star trackers on the spacecraft. The pointing 
for each bolometer is computed by combining the attitude with 
the location of the bolometer in the focal plane, as determined 
by planet observations (see below). Time-ordered informations 
of raw bolometer data are first processed to produce cleaned 
timelines and to set flags to mark data we do not currently fit. 
This TOI processing includes (1) signal demodulation and filter- 
ing, (2) deglitching, which flags the strong part of any glitch and 
subtracts the tails, (3) conversion from instrumental units (volts) 
to physical units (watts of absorbed power, after a correction 
for the weak non-linearity of the response), (4) decorrelation of 
thermal stage fluctuations, (5) removal of the systematic efl'ects 
induced by 4 K cooler mechanical vibrations, and (6) deconvo- 
lution of the bolometer time constant. Focal plane reconstruction 
and beam shape estimation are made using observations of 
Mars. The simplest description of the beams, an elliptical 
Gaussian, leads to full- width at half-maximum (FWHM) values, 
Os, given in Table 3 of the HEP {i.e., 7.08', 4.71', 4.50', 4.72' 
and 4.42' from 143 to 857 GHz, with an uncertainty between 
0.12'and 0.28'). From the cleaned TOI and the pointing, 
channel maps were computed using all the bolometers at a 
given frequency. The path from TOI to maps in the HFI DPC is 
schematically divided into three steps: ring-making, ring off'set 
estimation, and map-making. The first step combines the data 
within a stable pointing period, during which the same circle 
on the sky is scanned repeatedly to create rings with higher 
signal-to-noise ratio, taking full advantage of the redundancy 
of observations provided by the Planck scanning strategy. The 
low-frequency component of the noise is accounted for in a 
second step by using a destriping technique that models this 
component as an off'set of the ring values. Finally, cleaned 
maps are produced by coadding the off'set-corrected rings. The 
maps are produced in Galactic coordinates, using the HEALPix 
pixelisation scheme (see http://healpix.jpl.nasa.gov and Gorski 
et al. (2005)). Photometric calibration is performed either at ring 
level (using the CMB dipole) for the lower frequency channels 
or at the map level (using FIRAS data) for the higher frequency 
channels (545 and 857 GHz). The absolute gain calibration of 
the HFI Planck maps is known to better than 2% for the lower 
frequencies (143, 217 and 353 GHz) and 7% for the higher 
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Figure!. Wiener filter applied to the 143 GHz map for CMB subtrac- 
tion. The filter essentially cuts out high multipoles where the CMB-to- 
noise ratio of the 143 GHz map is low. Whereas this filter has to be 
known for estimating and subtracting the contribution of the residual 
CMB and 143 GHz instrument noise to the power spectrum of CMB- 
cleaned channels, the exact value of the filter is not really critical. 



frequencies (545 and 857 GHz), as summarised in the HEP 
Table 3. Inter-calibration accuracy between channels is better 
than absolute calibration. 

We made use of the so-called DX4 HFI data release, a dataset 
from which the CMB has not been removed. We used the 217, 
353, 545, and 857 GHz channels for the CIB analysis and the 
143 GHz channel for CMB removal. Maps are given either in 
MJysr"^ (with the photometric convention y/y =constant^) or 
yuKcMB , the conversion between the two was exactly computed 
using the bandpass filters (see Planck HFI Core Team 201 la). 

2.2. Extragalactic fields with high angular resolution Hi data 

Although Planck is an all- sky survey, we restricted our first CIB 
anisotropy measurements to a few fields at high Galactic lati- 
tude to minimize the Galactic dust contamination. The choice of 
the fields was driven by the availability of Hi data at an angular 
resolution close to HFI. 

The 21 -cm Hi spectra used here were obtained with the 100- 
m Green Bank Telescope (GBT) over the period 2005 to 2010. 
Details of this high-latitude survey are presented by Martin et al. 
(in prep). The total area mapped is about 825 deg^. 

The spectra were taken with on-the-fly mapping. The pri- 
mary beam of the GBT at 21 cm has a FWHM of 9. T, and the in- 
tegration time (4 s) and telescope scan rate were chosen to sam- 
ple every 3.5', more finely than the Nyquist interval, 3.86'. The 
beam is only slightly broadened to 9.4' in the in-scan direction. 
Scans were made moving the telescope in one direction (galactic 
longitude or right ascension), with steps of 3.5' in the orthogonal 
coordinate direction before the subsequent reverse scan. 

Data were recorded with the GBT spectrometer by in-band 
frequency switching, yielding spectra with a velocity coverage 
-450 < Vlsr < +355 kms"^ at a resolution of 0.80 kms"^ 
Spectra were calibrated, corrected for stray radiation, and placed 
on a brightness temperature (7^) scale as described in Blagrave 
et al. (2010), Boothroyd et al. (in prep), and Martin et al. (in 

^ The convention y/v=constant means that the MJy sr~^ are given for 
a source with a spectral energy distribution ly oc For a source with 
a diff'erent spectral energy distribution a colour correction has to be ap- 
plied (see Planck HFI Core Team 201 la). 
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prep). A third-order polynomial was fitted to the emission-free 
regions of the spectra to remove any residual instrumental base- 
line. The spectra were gridded on the natural GLS (SFL) projec- 
tion to produce a data cube. Some regions were mapped two or 
three times. With the broad spectral coverage, all Hi components 
from local gas to high- velocity clouds are accessible. 

We selected from this GBT cirrus survey the six faintest 
fields in terms of Hi column densities. Their main characteris- 
tics are given in Table 1 and the IRAS lOO/zm maps are shown 
in Fig. 1. They all have very low dust contamination and conse- 
quently Hi column densities, including the faintest all- sky sight 
line (referenced as LH2 in the Table). The field areas are between 
16 to 25deg^ for a total coverage of about 140 deg^. Going to 
higher average Hi column densities (N(Hi)> 2.5 x 10^^ cm"^) is 
not recommended because dust emission associated with molec- 
ular gas starts to contaminate the signal (see Fig. 9 of Planck 
Collaboration 201 It) and Hi is no longer a good tracer of dust 
emission. 

The HEALPix HFI maps were reprojected onto the small 
Hi maps by binning the original HEALPix data (sampled with 
HEALPix Nside of 2048, corresponding to a pixel size of 1.72') 
into Hi map pixels (pixel size 3.5' for all fields). An average of 
slightly more than four HEALPix pixels were averaged for each 
small map pixel. 
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2.3. Removing the bright sources from HFI maps 

We removed from the maps all sources listed in the Planck 
Early Release Compact Source Catalog (ERCSC) (Planck 
Collaboration 2011c). This represents only a few sources per 
field (if any), but the bright source removal is important for 
both power spectrum analysis and CMB map construction. It 
is also important to know the flux limit to compute the radio 
and dusty galaxy shot-noise contribution to the power spectra. 
Since our fields have roughly the same (very low) dust contami- 
nation, source detection is not limited by cirrus. Indeed, the flux 
cut is set by extragalactic source confusion at high frequencies 
and CMB contamination at low frequencies. The same flux cut 
can therefore be applied to all our fields. We took the minimum 
ERCSC flux densities in our fields as the flux cuts. They are 
given in Table 3. 

In practice, point source removal is performed in the original 
HFI HEALPix data prior to reprojection. For each source, a disc 
of size equal to the FWHM of the beam centred on the source 
position is blanked. Holes caused by missing data are then filled 
by a gap-filling process, which interpolates/extrapolates into the 
hole the values of neighbouring pixels. 
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2.4. Removing the CMB contamination from HFI maps 

Cosmic Microwave Background anisotropies contribute signifi- 
cantly to the total HFI map variance in all channels at frequen- 
cies up to and including 353 GHz. The detection and characteri- 
sation of CIB anisotropies at these frequencies requires the sep- 
aration of the contribution from the CMB. 

The present work focuses on very clean regions of the sky, 
for which Galactic foregrounds are very faint and are monitored 
using ancillary Hi observations. To remove the CMB in the 
fields retained for our analysis, we used a simple subtraction 
technique. While this simple method could be improved in 
future, it enables us to reliably evaluate CMB residuals, noise 
contamination, and to propagate errors due to imperfect instru- 
mental knowledge. It also guarantees that high-frequency CIB 



Figure 3. Power spectra of the different components for field SP (the 
figure is similar for the other fields). Power spectra of the 217, 353, 
545, and 857 GHz Planck maps (continuous black line) are compared 
to the noise power spectra (diamonds), to the CMB -cleaned power spec- 
tra (red), and to the CMB- and interstellar dust-cleaned power spectra 
(green). In this plot signal power spectra have not been corrected for 
the beam window function. Noise power spectra are computed using 
half-pointing period maps, as explained in Sect. 3.3. 



anisotropy signals will not leak into lower frequency, CMB -free 
maps. 

We removed the CMB contamination in the 217 and 
353 GHz channels by subtracting a CMB template obtained from 
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Figure 4. Maps of the 26 deg^ of the Nl field, from left to right: 217, 353, 545 and 857 GHz. From top to bottom: raw HFI maps; CMB- and ERCSC 
source-cleaned maps; residual maps (CMB-, sources-, and cirrus-cleaned); residual maps smoothed at 10' to highlight the CIB anisotropies. The 
joint structures clearly visible (bottom row) correspond to the anisotropies of CIB. Residual point sources are also visible. They have fluxes lower 
than the fluxes of the ERCSC removed sources. They have no impact on our analysis. 



the lower frequency data. We modelled the data for each fre- 
quency y as 



,CMB 



+ a 



,CIB 



(1) 



where x^miy) represents the channel data at frequency v, a^^^ 
the CMB map, ci^^iy) the CIB map and is a noise term 

comprising (if needed) any other astrophysical contaminant. The 
effect of the beam was accounted for with a (channel dependent) 
multiplicative factor, Z?^(y) (see Sect. 3.2). For the purpose of 
CMB removal, b^iy) was obtained from the Gaussian best-fit to 
the efifective HFI beam of the channel maps. 



At 100 and 143 GHz, we assumed that in the fields of inter- 
est only CMB and noise is present. Cosmic infrared background 
anisotropies are very small, and in the selected fields the contam- 
ination by other sources (e.g., cirrus) is negligible (see Sect. 3.4). 
In principle, both channels can be used to make a template of 
CMB emission. However, the 100 GHz channel is significantly 
less sensitive than 143 GHz and has an angular resolution two 
times worse than the 217 and 353 GHz channels. Therefore, we 
only used the 143 GHz channel as a CMB template. We cor- 
rected the 217 GHz maps for CMB contamination as follows 



6 



Planck Collaboration: CIB anisotropics with Planck 



= h{v2xi) [af^iyiii) + (1 - 

^Kvi43) 



(2) 



where W£ is a Wiener filter, designed to minimize the total 
contamination of the new map, yiyiii), by CMB and instrument 
noise. The 353 GHz map is corrected from CMB contamination 
in a similar way. Note that this cleaning was performed on a large 
region comprising all small fields used in the present analysis. 
The Wiener filter is obtained as 



We = 



(3) 



where C™^ is the current best-fit CMB model spectrum, and 
^^^(^143) is the power spectrum of the 143 GHz map. The Wiener 
filter Wi is close to 1 at low and close to at large i (see Fig. 
2). Note that errors on the beam estimate, on the assumed CMB 
power spectrum, or on the estimation of the 143 GHz power 
spectrum would result in sub-optimal filtering rather than in bi- 
ases. We checked that the CMB remaining in the CMB -cleaned 
maps does not change significantly with diff'erent assumptions 
leading to diff'erent w^. 

Errors in photometric calibration between channels are a 
problem. Although these errors are estimated to be small (2% 
at 143, 217, and 353 GHz), they may result in residual CMB at 
low i. They are accounted for in the processing, as detailed in 
Sect. 4.2.1. 

Fig. 3 shows the HFI power spectra of the raw and CMB- 
cleaned maps for one of the six fields. The CMB correction 
is very large at 217 GHz: the residual is a factor -100 below 
the raw power spectrum at ^ ^ 430 (it is a factor ~2 below at 
353 GHz). Note that whereas this illustrates the eff'ectiveness of 
CMB removal, it is also a source of worry about the impact of 
relative calibration errors for the 217 GHz channel. However, 
the power spectrum after CMB cleaning is -1% of the original 
map power spectrum only for € < 600. Cosmic microwave 
background-cleaned maps are shown on Fig 4. We see that the 
CMB has been efficiently removed. 

Finally we remark that an alternative method of removing 
CMB contamination, based on an internal linear combination of 
frequency maps and a needlet analysis (Delabrouille et al. 2009), 
was extensively studied and used in some of the Planck early 
papers, but it was not well suited to our purposes. The method 
tended to perform well over large patches of sky but left visible, 
large-scale residuals in the sky patches of interest, and had leak- 
age between the faint CIB and the CMB when other components 
(noise and Galactic cirrus) are present. 



Nl 


217 GHz 


353 GHz 


545 GHz 


857 GHz 


217 GHz 


\ 


0.56 


0.53 


0.49 


353 GHz 




1 


0.84 


0.77 


545 GHz 






1 


0.91 


Bootes 1 


217 GHz 


353 GHz 


545 GHz 


857 GHz 


217 GHz 


1 


0.44 


0.39 


0.39 


353 GHz 




1 


0.75 


0.74 


545 GHz 






1 


0.89 



Table 2. Pearson correlation coefficient between CIB anisotropy maps 
(values are given for the Nl and Bootes 1 fields to illustrate the range 
of coefficients). The high-frequency maps are highly correlated. A 
decorrelation is seen when going to lower frequencies. We interpret this 
decorrelation as reflecting the redshift distribution of CIB anisotropies 
(see text, Fernandez-Conde et al. (2008) and Penin et al. (201 la)). 



Hi components -The Hi data in each field show different veloc- 
ity components: a local component, typical of high-latitude Hi 
emission, intermediate- velocity clouds (IVCs) and high- velocity 
clouds (HVCs). These clouds are typically defined as concentra- 
tions of neutral hydrogen at velocities inconsistent with a simple 
model of differential Galactic rotation. The distinction between 
IVCs and HVCs is loosely based on the observed radial veloc- 
ities of the clouds; IVCs have radial velocities with respect to 
the local standard of rest (LSR) of 30 < IVlsrI < 90km s-\ 
while HVCs have velocities IVlsrI > 90 km s"^. High- velocity 
clouds might be infalling clouds fueling the Galaxy with low- 
metallicity gas, whereas IVCs might have a Galactic origin (e.g. 
Richter et al. 2001). For each ffeld, we constructed integrated Hi 
emission maps of the three velocity components. The selection 
of the velocity range for each component was based on inspec- 
tion of the median 21 -cm spectrum and of the rms 21 -cm spec- 
trum (i.e., the standard deviation of each channel map). It is fully 
described in the Planck Collaboration (201 It). The Hi maps were 
then converted to Hi column density using the optically thin ap- 
proximation: 



7V(HI)(jc,3;) = 1.823 x 10^^ ^ T^,{x,y ,y)6y , 



(4) 



where Tb is the 21 -cm brightness temperature and v the ve- 
locity. Corrections have been applied for opacity (see Planck 
Collaboration 201 It), they are lower than 5% for our CIB ffelds. 
As illustrated in Fig. 5, the different ffelds have clearly distinct 
Hi contributions, with, e.g., no local component in the direction 
of the AG ffeld. 



Hi-dust correlation -To remove the cirrus contamination from 
the HFI maps, we need to determine the far-IR to mm emission 
of the diff'erent Hi components identiffed with the 21 cm obser- 
vations. We assumed that HFI maps, /y(x,y), at frequency v can 
be represented by the following model 



(5) 



2.5. Removing the cirrus contamination from HFI maps 

From 100 jim to 1 mm, at high Galactic latitude and outside 
molecular clouds a tight correlation is observed between far- 
infrared emission from dust and the 21 -cm emission from gas^ 
(e.g. Boulanger et al. 1996; Lagache et al. 1998). Hi can thus be 
used as a tracer of cirrus emission in our ffelds, and indeed it is 
the best tracer of diffuse interstellar dust emission. 



The Pearson correlation coefficient is > 0.9 (Lagache et al. 2000). 



where N^^^ix^y) is the column density of the ith Hi component, 
is the far-IR to mm - Hi correlation coeflficient of component 
/ at frequency v and Cy(x, y) is a residual. The correlation coef- 
ff cients (often called emissivities) were estimated using a 
minimization given the Hi and HFI data and the model (Eq. 5). 
Although the Hi column densities of the diff'erent velocity com- 
ponents are quite similar (see Fig. 5), the emissivities may vary 
by factors of more than 10 between the local/IVC and HVC (see 
Planck Collaboration 201 It), so it is important to consider them 
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N(HI) 
2.6 





MJy/sr 
1.29 



Figures. Hi and dust maps for two fields: SP (top) and AG (bottom). The first two maps on the left show the Hi components (Local and IVC for 
SP, IVC and HVC for AG), the third maps show the 857 GHz emission associated with Hi (2^ Q^yA^Hj) maps on the right side show the HFI 

857 GHz maps. Those HFI maps have been convolved by the GBT beam to allow a better comparison by eye. Hi maps are given in units of 10^*^ 
atoms cm~^. Note the correlation of the dust emission with the diff'erent Hi velocity components and its variation from field to field. 



separately. The emissivities can be used to characterise the opac- 
ity and temperature of the dust emission in the diflTerent compo- 
nents. This is beyond the scope of this paper, but it is extensively 
discussed in the Planck Collaboration (201 It). 

Cirrus contamination removal -We removed from the HFI 
maps the Hi velocity maps multiplied by the correlation coeffi- 
cients. Maps are shown in the last two columns of Fig. 5 for two 
fields. The removal was made at the HFI angular resolution, even 
though the Hi map is of lower resolution. This is not a problem 
because cirrus, with a power-law power spectrum (Miville- 
Deschenes et al. 2007), has negligible power between the GBT 
and HFI angular resolutions, in comparison to the power in the 
CIB. 



Residual maps and power spectra -Fig. 3 shows the HFI 
power spectra before and after the dust removal in the SP field. 
Cirrus removal has more impact for the two high-frequency 
channels. At 217 GHz, the correction is very small (13% at 
^=500). This method of using Hi data to remove the cirrus con- 
tamination from power spectra has also been successfully ap- 
plied by Penin et al. (2011b) at higher frequencies than ours, 
where the cirrus contamination is higher. The authors have been 
able to isolate precisely the CIB anisotropies power spectra at 
1875 and 3000 GHz with Spitzer and IRAS/IRIS, in the Nl field. 

The residual maps at the HFI angular resolution are shown 
in Fig. 4 for the Nl field. We clearly see that the cirrus has 
been efficiently removed. The bottom row shows the residual 
maps, smoothed at 10'. Common structures, corresponding to 
the CIB anisotropies, are clearly visible at the four frequencies. 
Table 2 gives the Pearson correlation coefficients between the 
CIB anisotropy maps. They are about 0.9 between the 545 and 
857 GHz maps and 0.5 between the 217 and 857 GHz CIB maps. 




V [GHz] 

Figure 6. Contribution to the CIB per redshift slice, extracted from 
Bethermin et al. (2011). The black solid line is the CIB spectrum pre- 
dicted by the model. The contribution to the CIB from < z < 0.3, 
0.3<z<l,l<z<2 and z > 2 galaxies is given by the red short- 
dashed, green dot-dashed, blue three dot-dashed and purple long-dashed 
lines, respectively. Lower limits coming from the stacking analysis at 
100 yL/m, 160 pm (Berta et al. 2010), 250 pm, 350 pm, 500 pm (Marsden 
et al. 2009), 850 pm (Greve et al. 2010) and 1.1 mm (Scott et al. 2010) 
are shown as black arrows. The black diamonds give the Matsuura 
et al. (2011) absolute measurements Wiih AKARI. The black square the 
Lagache et al. (2000) absolute measurements with DIRBE/WHAM and 
the cyan line the Lagache et al. (2000) FIR AS measurement. 



The decrease when the frequency diflTerence between the maps 
is larger is expected because the contribution of high-redshift 
galaxies to the CIB (and its anisotropies) increases with wave- 
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Figure?. A number of recent models of dusty-galaxy evolution and 
their associated shot noise for different flux cuts at 857 GHz. Top: 
Comparison of the models with the Herschel and BLAST differen- 
tial numbers counts. Models are from Lagache et al. (2004); Negrello 
et al. (2007); Le Borgne et al. (2009); Patanchon et al. (2009); Pearson 
& Khan (2009); Vahante et al. (2009); Bethermin et al. (2011); 
Franceschini et al. (2010); Lacey et al. (2010); Marsden et al. (2011); 
Rowan-Robinson (2009); Wilman et al. (2010). Data points are from 
Oliver et al. (2010); Bethermin et al. (2010b); Glenn et al. (2010). 
Bottom: Shot-noise level as a function of the flux cut for the same mod- 
els (same colour and line coding between the two figures). The verti- 
cal and horizontal continuous dark lines show the Planck flux cut and 
shot-noise level from Table 3, respectively. The Bethermin et al. (201 1) 
model is shown by the continuous dark line. This figure shows that 
models predicting a very high shot noise (e.g. continuous and dashed 
light-blue, red-dashed, continuous and dashed dark-blue lines) are in- 
compatible with the measured number counts. 



length. This is illustrated in Fig. 6, extracted from Bethermin 
et al. (201 1), where we show the redshift distribution of the CIB. 
The redshift distribution of correlated CIB anisotropies is dis- 
cussed in Fernandez-Conde et al. (2008, 2010) and Penin et al. 
(2011a). 



3. Astrophysical and instrumental components of 
residual HFI maps power spectra 

Once the CMB and cirrus have been removed, there are three 
main astrophysical contributors to the power spectrum at the HFI 
frequencies: two from dusty star-forming galaxies (with both 
shot noise, C^''^°^(y), and clustering, cY^''''\v), components), 
and one from radio galaxies (with only a shot-noise component. 



° (y)^ the clustering of radio sources being negligible, see 
Hall et al. (2010)). If the instrument noise and the signal are not 
correlated, the measured power spectrum Q(y) is 



Q(y) = /72(y)[cf^-^(y) + Cf^^\y) + C^'^^^^(y)] 



(6) 



where b^iy) is the beam window function, and N^iy) the power 
spectrum of the instrument noise. Note that we neglect here the 
Sunyaev-Zeldovich (SZ; Sunyaev & Zeldovich 1980) contribu- 
tion to the power spectra. Extrapolation of SPT (Lueker et al. 
2010) and ACT (Dunkley et al. 2011) constraints show that SZ 
is negligible compared to CIB anisotropies at y > 217 GHz. Our 
goal is to accurately measure C^'^^^^^(y), which we present and 
extensively discuss in Sect. 5. We begin by discussing all other 
components of Eq. 6 in this section. 

3.1. Shot noise 

The shot noise arises from sampling of a background composed 
of a finite number of sources. We assumed the distribution is 
Poisson, so that its power spectrum is independent of L If we 
identify and remove all sources brighter than 5'cut, the shot noise 
from the remaining sources fainter than 5'cut is given by (e.g., 
Scott & White 1999) 



/^shot 



f 

Jo 



dS 



dS, 



(7) 



where S is the source flux and dN/dS the differential number 
counts. These counts can be directly measured or derived from 
evolution models of the relevant population of galaxies (dusty, 
star-forming and radio galaxies in our case). 



3.1 .1 . Star-forming, dusty galaxy shot noise, C^''^°^ 

We used the recent model of Bethermin et al. (2011) to com- 
pute the IR galaxy shot-noise power. This is an updated version 
of the Lagache et al. (2004) model that better reproduces new 
observational constraints (e.g., from Herschel). This new, em- 
pirical model uses the same galaxy spectral energy distribution 
(SED) templates as Lagache et al. (2004), but a fully paramet- 
ric evolution of the luminosity function. The parameters of the 
model were determined by fitting the infrared/sub-mm number 
counts, and some mid-IR luminosity functions, with a Monte- 
Carlo Markov Chain (MCMC). More details on the model are 
given in Appendix B. The derived shot-noise power is given 
in Table 3, with uncertainties computed from the MCMC. The 
quoted numbers include statistical and photometric calibration 
uncertainties. This model has less energy output at high redshift 
(z - 2) and consequently lower shot-noise power at long wave- 
length than the Lagache et al. (2004) model. The shot noise lev- 
els depend on the flux cut, which itself has an uncertainty linked 
to the flux uncertainty in the ERCSC. If we change the flux cut 
^cut by 30% in Eq. 7 based on the uncertainty in ERCSC fluxes, 
the power spectra change by less than 5% at all frequencies (and 
less than 1% at 217 GHz). 

As we will discuss in Sect. 5, the dusty galaxy shot- 
noise level will be a major factor in the interpretation of CIB 
anisotropy power spectra. Because we are obtaining this value 
from a model, not measuring it directly in this paper (see 
Sect. 5), we briefly discuss here the constraints on the model and 
the plausible range of values using the 857 GHz channel as an 
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example (the same conclusions are reached for the other Planck 
channels). Fig. 7 shows a compilation of models from the lit- 
erature superimposed on the latest number counts observed by 
BLAST and Herschel and the expected shot noise as a function 
of 5'cut- First we see, as stated above, that a small variation in 5'cut 
leads to only a small variation in shot-noise power. Second, we 
see that the highest shot-noise level is around 13,500 Jy^ sr"\ 
a factor ~ 2.3 above our nominal value, but it comes from a 
model that overestimates the observed number counts by a large 
factor (3-4 for 50 < 5 < 300 mJy). Models that agree rea- 
sonably well with the number counts have a shot-noise level 
below SOOOJy^sr-^ The Bethermin et al. (2011) model has 
the lowest shot noise. However, it is currently the model that 
best reproduces all available constraints, from the mid-infrared 
to the millimeter, including the differential contribution of the 
5" 24 > 80 yu Jy sources to the CIB as a function of redshift, which 
is a difficult observation to predict. Eventually, the shot noise de- 
rived from this model agrees well also with the HerscheljSV^RE 
measurements in the Lockman- SWIRE field, when none of the 
point sources is removed (Amblard, private communication), as 
detailed in Sect. 5.3. 



1.0 




1000 1200 1400 1600 1800 2000 2200 2400 
Multipole I 



Figure 8. Effective beam window functions (bi) from FICSBell (black) 
and FEBeCoP (red) at 545 GHz (see Sect. 3.2 for more details). The six 
FEBeCoP beam window functions from each field are superimposed 
(red lines). Also shown for comparison is the Gaussian beam with a 
FWHM of 4.72' + 0.2' (green lines), which is the equivalent FWHM of 
the beam determined on Mars. 



3.1 .2. Radio galaxy shot noise, Cf 

The shot-noise power from radio galaxies is subdominant to 
that from dusty sources at the frequencies relevant to CIB 
anisotropy analysis. The radio galaxy shot-noise power can be 
estimated from the model of de Zotti et al. (2005). At frequencies 
< 100 GHz, the model agrees with the source counts computed 
using the extragalactic radio sources from the ERCSC. At 143 
and 217 GHz, and for fluxes below 300 mJy {i.e., the case listed 
in Table 3) the de Zotti et al. (2005) model agrees with the source 
counts of Vieira et al. (2010). At higher fluxes the model needs 
to be scaled to reproduce the number counts obtained using the 
ERCSC. The estimated scaling factors are 2.03 and 2.65 at 143 
and 217 GHz, respectively (see Planck Collaboration 20111). At 
even higher frequencies the number counts by the ERCSC are 
no longer complete. We therefore use the 217 GHz scaling factor 
to set upper limits for the shot noise. It is negligible compared 
to cf^^^^ at these frequencies (see Table 3). Changing the flux 
cut by 30% affects the shot noise by 30%, but because the radio 
contribution is subdominant at the frequencies relevant for CIB 
anisotropy analysis, it has little impact on our results. 

3.2. The beam window function, b^iv) 

Because the HEI beams are not azimuthally symmetric, the scan- 
ning strategy has to be taken into account in modelling the effec- 
tive beam response. We used two different methods to compute 
the effective beam: EEBeCoP and EICSBell. With EEBeCoP, we 
computed one effective beam per field, with EICSBell, one ef- 
fective beam for the entire sky. 

FICSBell -The EICSBell method (Hivon et al, in prep) gener- 
alizes the approach of Hinshaw et al. (2007) and Smith et al. 
(2007) to polarization and to include other sources of systemat- 
ics. The different steps of the method used for this study can be 
summarized as follows: 

1. The scanning-related information (i.e., statistics of the orien- 
tation of each detector within each pixel) is computed first. 



and only once for a given observation campaign. The hit 
moments are only computed up to degree 4, for reasons de- 
scribed below. 

2. The (Mars-based) beam map or beam model of each detec- 
tor, d, is decomposed into its spherical harmonic coefficients 

bi = J dr BAr)Y,s(r), (8) 

where B^ir) is the beam map centred on the North pole, 
and Y{s(y) is a spherical harmonic. Higher s indices de- 
scribes higher degrees of departure from azimuthal symme- 
try and, for HEI beams, the coefficients b^^ are decreasing 
functions of s at most £ considered. It also appears that for 
£ < 3000, the coefficients with |^| > 4 account for < 1% 
of the beam throughput. For this reason, only modes with 
1^1 < 4 are considered in the present analysis (Armitage- 
Caplan & Wandelt (2009) reached a similar conclusion in 
their analysis of Planck-LFl beams). 

3. The b^^ coefficients computed above are used to generate ^- 
spin weighted maps for a given CMB sky realisation. 

4. The spin- weighted maps and hit moments of the same or- 
der, s, are combined for all detectors involved, to provide an 
"observed" map. 

5. The power spectrum of this map can then be computed, and 
compared to the input CMB power spectrum to estimate the 
effective beam window function over the whole sky, or over 
a given region of the sky. 

Monte-Carlo (MC) simulations in which the sky realisations are 
changed can be performed by repeating steps 3, 4, and 5. The 
impact of beam model uncertainties can be studied by including 
step 2 into the MC simulations. 

FEBeCoP -As mentioned in Sect. 2.1, map making reduces 
time-ordered data to pixelised maps. Each pixel of a map rep- 
resents a convolution of the true sky with the combined effect 
of scanning beam and scan pattern. EEBeCoP computes this 
combination of beams and scans — the effective beams — as is, 
in the pixel space. The EEBeCoP methodology and algorithm 
has been described in Mitra et al. (2011), and Planck HEI Core 
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Frequency (GHz) 


143 


217 


353 


545 
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Flux cut (mJy) 


245 


160 


325 


540 
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IR shot noise ^ 


1.4+0.3 


12.2+2.9 


138+22 


1150+92 


5923 +367 


(Jy2 sr-i) 












Radio shot noise^ 


7.1 


4.0 


<3.4 


<5.7 


<7.4 


(Jy2 sr-i) 












IR shot noise ^ 


(1.0 + 0.2) X 10"^ 


(5.3 + 1.2) X 10-^ 


(1.7 + 0.3) X 10-^ 


0.34 + 0.03 


1187 + 74 














Radio shot noise^ 


5.2 X 10"^ 


1.7 X 10"^ 


<4.1 X 10"^ 


< 1.7 X 10"-^ 


< 1.5 















Table 3. Flux cut from the ERCSC for our six fields, and the shot-noise power for dusty and radio galaxies appropriate to those cuts (see text). 
Values for shot noise^ are derived from the dusty galaxy evolution model of Bethermin et al. (201 1), while those for shot noise^ are from the radio 
galaxy evolution model of de Zotti et al. (2005) (see text for more details). 
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Figure 9. Three independent noise-power- spectrum measurements in 
the SP field at 353 GHz: red continuous line, half pointing period; green 
dashed, surveys I and II; black dot-dashed, half focal plane array). 



Team (2011a). Below, we list for completeness the essential 
steps made in computation of the beam window functions: 

1 . For each pixel / in the map (or CIB field) we computed the 
Fourier-Legendre transform, B^, of the pixel space eflTective 
beams Bi(Ci) using the formula 

= r dClP,iCli-Ci)Bi(Ci), (9) 

where O/ is the direction vector of the centre of the /th pixel 
on the sky, represents Legendre polynomials of order I 
and the integration is performed over the (small) solid angle 
AO/, outside which the beam can be taken as zero. This for- 
mula can be readily transformed to a discretised form with a 
careful correction for the "pixel window function" as 

j 

where the summation is over pixels that fall inside the beam 
solid angle AO/, Dpix is the area of each (equal area) pixel 
and is the pixel window function that compensates for 
the systematic error that is introduced when integration over 
a pixel is replaced by the value of the integrand at the pixel 
centre times the area of the pixel. 

2. We then computed at uniformly sampled directions in 
each field to find the average window functions. The samples 



were chosen as the HEALPix pixel centres at a coarser res- 
olution (A^side = 128) to ensure uniform sampling. Thus we 
obtained the average window functions for each frequency 
and field. 

3. To validate the average window functions obtained using 
the above prescription, we performed Monte-Carlo simula- 
tions separately for each field and each frequency. We sim- 
ulated 16 realisations of the sky starting from a oc angu- 
lar power spectrum, which are convolved in two ways - (1) 
with FEBeCoP-generated eflTective beams in pixel space and 
(2) with analytical Gaussian beam in harmonic space for a 
beam size appropriate for the given frequency channel. The 
convolved maps were then "masked" using a function that is 
unity in the given field and smoothly (in ~ 25% of field ra- 
dius) goes to zero outside the field. Finally, we computed the 
ratio of the angular power spectra of these two maps, mul- 
tiplied the ratio by the theoretical window function for the 
same beam size and averaged over the realisations. Though 
these "transfer functions" suflTer from ringing eflfects often 
seen in Fourier transforms of a narrow function, they wiggle 
around the average window functions, confirming the valid- 
ity of the latter. 

Fig. 8 shows the FICSBell and FEBeCoP eflfective beams at 
545 GHz. Also shown is the Gaussian beam with a FWHM of 
4.72' ± 0.2r. This is the average FWHM of the scanning beam, 
determined on Mars obtained by unweighted averaging the indi- 
vidual detectors FWHM. Each FWHM is that of the Gaussian 
beam, which would have the same solid angle as that deter- 
mined by using a full Gauss-Hermite expansion on destriped 
data (see Planck HFI Core Team 2011a, for more details). We 
see a quite good agreement between the FICSBell, all-sky and 
FEBeCoP, small-field eflfective window functions, with a 2% dif- 
ference at ^ ^ 2000, the highest I that will be considered for CIB 
anisotropy analysis (see Sect. 4). We also see from the figure that 
the error on the input scanning beam is larger than this diflference 
and will dominate the uncertainties at high i and high frequency 
(Sect. 4.2.2). Below we will use the FEBeCoP window functions 
because they are exactly computed for each of our fields. 

3.3. Instrument noise, Nf(y) 

We can use three diflTerent jack-knife diflTerence maps to derive 
noise power spectra: maps made from the first and second halves 
of each pointing period (a half-pointing period is of the order of 
20 minutes), maps made using half of the focal plane array, and 
maps using the two diflTerent surveys (surveys I and II). In each 
case the noise power spectrum, A^^, is obtained by measuring the 
power spectrum of the diflTerence maps. The three methods give 
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1, 3, and 4% at 217, 353, 545 and 857 GHz, respectively. Fig. 10 
shows the noise power spectra for all fields. They are nearly flat, 
the deviation from flatness is caused by the eff'ect of deconvolu- 
tion from the instrumental response at high frequency and resid- 
ual low-frequency noise. Removing the ERCSC sources has no 
impact on the noise determination. 

Fig. 3 shows the noise power spectra compared to the HFI 
map power spectra for one illustrative field. We see that we have 
a very high signal-to-noise ratio. At 545 and 857 GHz, the signal 
is dominating even at the highest spatial frequencies. At 217 and 
353 GHz, the residual signal (i.e., CMB- and cirrus-cleaned) is 
comparable to the noise at high f (f > 2000 - 2500 depending 
on the field). 

3.4. Additional corrections 

Two additional corrections linked to the CMB cleaning were 
made for the power spectra. First we removed the extra instru- 
ment noise that has been introduced by CMB removal: 



(11) 



with y equal to 217 or 353 GHz. A^^Cvhs) is the noise power spec- 
trum of the 143 GHz map. It is computed as the noise in the other 
frequency channels, using the half-pointing period maps, follow- 
ing Sect. 3.3. 

Second, owing to the lower angular resolution of the 
143 GHz channel compared to the 217 and 353 GHz, we also 
had to remove the CMB contribution that is left close to the an- 
gular resolution of the 217 and 353 GHz channels: 

^CMBres^^^ = C™^(y) X X b]{v) X [1 - W^f , (12) 

with Fp the pixel and reprojection transfer function (detailed in 
Sect. 4.1). 
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Figure 10. Instrument noise power spectra of the six fields obtained 
using half-pointing period maps. From top to bottom: 217, 353, 545 and 
857 GHz (continuous: NX, dotted: AG, dashed: SP, dash-dotted: Bootes 
1, long-dash: Bootes 2, dash-3 dotted: LH2). 



similar N^, as is illustrated for one frequency and one field in 
Fig. 9. We chose, however, to use the half-pointing period maps 
because (1) the two survey maps are only fully covered for the 
LH2 and SP fields and (2) there are only three bolometers at 545 
and 857 GHz, making half-focal plane maps less accurate. We 
also computed the noise power spectrum from the diflFerence be- 
tween the auto- and cross- power spectrum of the two half maps. 
In the range of interest, 1500 < £ < 2100, where the contribution 
from the noise becomes important, they agree at better than 0.5, 



Finally, we had to assess the level of the astrophysical 
components that were removed from or added to the 217 and 
353 GHz channels, using the filtered 143 GHz channel as a CMB 
template. Cirrus emission is highly correlated between 143, 217 
and 353 GHz channels. Consequently, filtered cirrus emission 
was removed from the 217 and 353 GHz. This has no impact 
on our CIB anisotropy analysis because this extra cirrus re- 
moval only modifies the emissivities, with no consequence on 
our residual maps (it should be understood for a further interpre- 
tation of the Hi-correlated dust emission, which is not the goal of 
this paper). We expect the shot-noise powers to be quite decor- 
related for the (143, 217) and (143, 353) sets of maps because 
the 143 GHz shot noise is dominated by radio sources, whereas 
the 217 and 353 GHz shot noise is dominated by dusty galax- 
ies (see Table 3). To have an idea of the maximal effect {i.e., 
perfect decorrelation between shot noise at 143, and 217, and 
353 GHz) we computed the contamination by the 143 GHz shot 
noise, sunmiing the contribution of the radio and dusty galaxies 
and following 



Q(v) = C^^°^(yi43) X 



Kv) 



^Kvi43) 



(13) 



The last term accounts for the filtering and "re-beaming" of the 
143 GHz map. The contamination is the highest in the 217 GHz 
channel. It is a factor 1 .2 and 120 smaller than the sum of the pre- 
dicted radio and dusty galaxies shot-noise powers at 217 GHz at 
i ^ 200, and 2000, respectively, but is equivalent at i ^1000. 
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Anyway, it is smaller by factors of 20, 2.9 and 325 than the 
CIB anisotropics at 217 GHz, at ^= 200, 1000, and 2000, respec- 
tively. Because this is the maximal contamination and because it 
is quite low (and completely negligible at high i), we did not 
apply any correction to the CIB anisotropy power spectra. 

We still have to consider the case of CIB -correlated 
anisotropics at 143 GHz. They have been marginally constrained 
at 150 GHz by SPT and ACT at high £. The power is < 5.2 x 
10" and < 9.8 x 10" at ^ = 3000 in Dunkley et al. 
(2011) and Hall et al. (2010), respectively. This contribution is 
also completely negligible compared to the signal at 217 GHz. 

In conclusion, we can ignore the CIB and cirrus components 
that are left in the CMB maps. 

4. Angular power spectrum estimation 

The angular power spectrum estimator used in this work is 
POKER (Ponthieu et al. 2011), which is an adaptation to the flat 
sky of the pseudo-spectrum technique developed for CMB anal- 
ysis (see e.g. MASTER, Hivon et al. 2002). In brief, POKER com- 
putes the angular power spectrum of the masked data (a.k.a., 
the pseudo-power spectrum) and deconvolves it from the power 
spectrum of the mask to obtain an unbiased estimate of the 
binned signal angular power spectrum. We summarize the main 
features of POKER in the following section and then detail how 
it was used to produce the final estimate of the CIB anisotropy 
power spectrum and its associated error bars. 

In the following, the power spectrum associated to CIB 
anisotropics will be denoted Q and its unbiased estimator in 
the flat-sky approximation Pif). As already suggested, this final 
estimate makes use of the power spectrum of the masked data. 
This so-called pseudo-power spectrum will be denoted P{i). In 
the flat-sky approximation, the standard angular frequencies la- 
belled by their zenithal and azimuthal numbers, usually called 
i and m respectively, are replaced by an "angular" wave-vector 
£\ its norm £ is equal to the zenithal number (see e.g., the ap- 
pendix of White et al. 1999). Finally, we will assume that CIB 
anisotropics arise from a statistically isotropic process. As is the 
case for the CMB, the CIB fluctuations are viewed as isotropic 
and homogeneous stochastic variables on the celestial sphere, 
leading to 

{a(i)a*(r)) = (27TfC(£)S^(e - r), (14) 

with a(i) the Fourier coefficients of CIB anisotropics. This as- 
sumption is theoretically reasonable, moreover, we checked that 
|a(/')p computed from our CIB maps does not depend on the di- 
rection of i. 

4.1. POKER 

The POKER implementation of the pseudo-spectrum approach 
uses the discrete Fourier transform (hereafter DFT). For a map 
of scalar quantity Djk (j, k denote pixel indices), it is defined as 

D,nn = Dj, X ,-2-0>"/^...«/^,)^ (15) 

Djk = YjDmnX e+^-X^/'v.+^^W), (16) 

m,n 

where Dmn is the set of discrete Fourier coefficients of 
Djk. For a given wave- vector i, labelled by the m and 
n indices, its corresponding norm is denoted by ^mn = 



UE + 00 




Figure 11. Contribution of residuals to the final CIB anisotropy esti- 
mate (illustrated here with the SP field at 353 GHz). The bias induced 
by each dust and CMB component is negligible compared to both the 
CIB anisotropy signal (black dots) and the statistical noise (in green, 
including cosmic variance on the noise estimate itself) 



(2n/Ae) ^(m'lN^f + in' INyf with m' = m (respectively n' and 
n) if m < Nx/2 and = Nx - m if m > Nx/2. The power spec- 
trum of the map is defined as the square-modulus of its Fourier 
coefficients, i.e., P(imn) = - 

The direct DFT of the masked data relates the true Fourier 
coefficients to the pseudo-V(^^xnQX coefficients of the signal 

m'n' 

in which W^'^, is a convolution kernel that depends only on the 

mask DFT coefficients. Replacing Dmn by Dmn in the definition 
of the power spectrum of a given map leads to the power spec- 
trum of the masked data (a.k.a. the pseudo-power spectrum). 
For a signal T plus noise N map, the ensemble averaged of the 
pseudo-spectrum tracing a statistically isotropic process, reads 

C(C'«') + <Wm«)>, (18) 

m'n' 

where we have introduced the total transfer function F^'n' ac- 
counting for the beam, the 'map-making' pixelisation eff'ects and 
reprojection from curved, HEALPix maps to flat, square maps. 
The beam transfer function is given by the beam power spectrum 
described in Sect. 3.2. The 'map-making' pixelisation efl'ects are 
described by the power spectrum of the pixel window function 
for full-sky maps provided by the HEALPix package (the ini- 
tial HEALPix maps are built with A^side = 2048 corresponding 
to a pixel size of 1.7'). As explained in Planck HFI Core Team 
(2011a), time-domain filtering is included as part of the scan- 
ning beam, such that any time-domain filtering eff'ects end up in 
the estimate of the beam instead of as part of F^'n' • Finally, each 
curved map with a 1.7' resolution is reprojected onto its tangent, 
flat space with a pixel size of 3.5'. This induces first a repixelisa- 
tion eff'ect because the output map is less resolved than the input 
one and second, a slight displacement of the pixel centres. The 
cumulative impact of 'image deformation' and 'repixelisation' 
is estimated via Monte-Carlo: we first generated a set of full- 
sky maps and computed the MC average of their pseudo-spectra. 
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This set of maps was then re-projected onto flat maps for which 
MC average of their pseudo- spectra in the flat-sky approxima- 
tion were computed. The ratio of the flat-sky pseudo- spectrum 
divided by the full-sky pseudo-spectrum gives a measurement of 
re-projection eff'ect. Note that those simulations have been per- 
formed assuming diff'erent shapes for the input angular power 
spectra. The derived reprojection transfer functions agreed per- 
fectly, which underlines the robustness of the approach. 

An unbiased estimate of Q is obtained by first subtract- 
ing the noise contribution and then deconvolving the mask and 

beam eflfects encoded in the convolution kernel \W^\\ F^'n'- 

I m,m I " 

For the sky coverage of the considered fields, the rapid oscil- 
lations of the convolution kernels introduce strong correlations 
between spatial frequencies and make its inversion numerically 
intractable. (Pseudo-)Power spectra are therefore estimated in 
frequency bands (labelled b hereafter). The binning operator is 



A. Given the mask, W, associated to our CIB map and the 
transfer function, Fm,n, compute and invert M^^/ as given by 
Eq. 24. The diff'erent fields (with a size at most 5.1° x 5.1°) 
are systematically embedded in a 10° x 10° square map. We 
use binary masks, i.e.,W = I for observed pixels (i.e., those 
kept in the analysis), and W = for unobserved pixels. 
The estimated power spectra are binned with a bandwidth 
of Ab = 200 and the first bin starting at ^ = 80"^. 

B. Derive the noise bias (Nt)^ given by first the instrument noise 
described in Sect. 3.3 and second the additional corrections 
given in Sect. 3.4. 

C. Compute the final estimate of Cf, from Eq. 25 and Cb = (Pb)- 

Uncertainties on Pb come from sampling variance, noise 
variance, astrophysical contaminants, and systematic eff'ects. In 
the following section, we present how we estimated each of these 
contributions from Monte-Carlo simulations. 



nmn 



otherwise 



fb+\ 
How 



(19) 4_2. Error bar estimation 



where Ab is the number of wave vectors tmn that fall into the bin 
b. The reciprocal operator that relates the theoretical value of the 
one-dimensional binned power spectrum Pb to its value at irnn is 



1 



low 



(20) 



otherwise 



For optimal results, the spectral index ft should be chosen to get 
£^C£ as flat as possible. For the CMB,j3 ^ 2 is the equivalent of 
the standard ^(^-F 1) pref actor. For the CIB anisotropics Q scales 
roughly as and we therefore adopted a binning with J3 = I. 
Nevertheless, we checked that our results were robust against the 
choice of JS: we simulated a power spectrum scaling as but 
reconstructed it assuming = 0, I and 2 in POKER. For each 
choice of fi, the estimated power spectrum perfectly agreed with 
the input one (for a more complete discussion, see Ponthieu et al. 
2011). 

The binned pseudo-power spectra is 



(21) 



m,neb 



and the CIB power spectrum is related to its binned value, Cb, 
via 

C(W)-2eL'C^' (22) 
With these binned quantities, Eq. 18 can be re-written as 

{Pb) = Yj^bb'Cb^^{Nb), (23) 



with 



M,,.= 2 Z ^ricr.'r^p-vGL'- (24) 



m,neb m',n'eb' 



An unbiased estimate of the binned angular power spectrum of 
the signal is thus given by 



Pb = Yj^bb'{Pb'-(^b'))^ 



(25) 



It is easily checked that (Pb) = C^. 

The complete recovery of the CIB anisotropy angular power 
spectra is therefore made in three steps: 



4.2.1. Statistical uncertainties 

As presented in Sect. 4.1, the uncertainties on Pb come from 
signal sampling variance, noise, and uncertainties on the sub- 
traction of CMB and Galactic dust. The first two are described 
by stochastic processes with known power spectra, the last two 
come from uncertainties in the weights applied to templates in 
the subtraction process at the map level. 

We developed the tools necessary to simulate maps given 
any input angular power spectrum for each field and frequency. 
The simulation pipeline consists of simulating maps given an 
input power spectrum (in the case of CIB, noise and CMB resid- 
ual) and maps of template residuals (conservative Gaussian ran- 
dom fractions of the templates). These maps are then combined 
and analysed by the power spectrum estimator. Each realisation 
provides an estimated power spectrum with the same statistical 
properties as our estimate on the data, and alltogether these sim- 
ulations provide the uncertainties on our estimate. The covari- 
ance matrix of Pb is 



Cbb' = ((Pb - (Pb)Mc) (Pb' - <^^'>Mc))mC ' 



(26) 



with (•)mc standing for Monte-Carlo averaging. The error bar on 
each Pb is 

o-pt = ^[^b. (27) 

and the bin-bin correlation matrix is given by its standard defini- 
tion 

= -J=- (28) 

^J^bb^b'b' 

Simulation pipeline -The simulated maps are 10° x 10°. They 
contain six components, accounting for the different ingredients 
supposedly present in the actual data maps: 

1. A CIB anisotropy component obtained from a random, 
Gaussian realisation of the CIB anisotropy power spectrum. 
As a model for such a spectrum, we used a fit to Pcm, es- 
timated from the data additionally multiplied by the power 
spectrum of the beam, pixel, and reprojection transfer func- 
tion; 



^ = 80 corresponds to the inverse of the largest angular scale con- 
tained in the considered fields. 
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2. a residual CMB component derived from a random, 
Gaussian realisation of the power spectrum given in Eq. 12 
using the known Wiener filter, beam diff'erences between 
143 GHz and other channels, and the WMAP best fit CMB 
temperature power spectrum; 

3. the instrument noise as derived in Sect. 3.3. Since the noise 
is slightly coloured, we simulate it using a fit of its measured 
power spectrum; 

4. extra instrument noise incurred by CMB removal using 
Eq. 1 1 as a model of its power spectrum. 

In addition to those four ingredients standing for signal and noise 
(a CMB residual being viewed as an extra-source of 'noise' from 
the viewpoint of CIB), we added the two foreground templates 
that are removed, with some uncertainties: 

5. A CMB map with a Gaussian uncertainty distributed with 2% 
and 3% standard deviation at 217 and 353 GHz, respectively 
(the CMB is negligible at higher frequencies compared to 
CIB and dust). The 2% and 3% are justified in Sect. 4.3; 

6. an Hi map with a 5%, 10% and 10% standard deviation 
for its emissivity (local, IVC and HVC components respec- 
tively), consistent with both the inter-calibration errors (see 
Sect. 4.3) and the emissivity errors computed by the Planck 
Collaboration (201 It) using Monte Carlo simulations. 

The analysis pipeline -The analysis pipeline works in four 
steps: 

1. A first set of 1,000 MC simulations of CMB residual and 
noise only is performed to assess first, the pseudo-power 
spectrum of the instrument noise and CMB residual used to 
debias the simulated data pseudo-spectrum and, second, the 



noise variance given by cr^^ = -^C 



inoise , 

bb • 



2. a second set of 1,000 MC simulations, including all the com- 
ponents, is performed. For a given simulated map, CMB and 
dust templates are removed assuming the estimated emissiv- 
ities of Sect. 2.5 for dust; 

3. the POKER algorithm is then applied to these 'foreground- 
cleaned' maps to obtain a final estimate of the CIB angu- 
lar power spectrum. In this step, the bias involved in Eq. 25 
contains the pseudo-power spectrum of the instrument noise 
model and of the CMB residual model as described in the 
simulation pipeline; 

4. the total error bars and bin-bin correlation matrix on Pi, are 
obtained as the RMS of 1,000 Monte-Carlo realisations of 
the simulation pipeline, as described in the previous para- 
graph and using Eqs. 27 and 28. 

The statistical uncertainties contain: 

A. sampling variance from CIB anisotropics and residual CMB 
as modelled in Eq. 12; 

B. noise variance from instrument noise and extra noise given 
by Eq 11; 

C. uncertainties on the CMB and Hi template subtraction. 

In this set of simulations and analysis, we assumed the beam 
profile as described in Sect. 3.2 and ignored potential beam un- 
certainties (see the next section for a discussion of this system- 
atic eff'ect). Below, we present the results obtained using the 
FEBeCoP-derived beam profiles, but the estimated power spec- 
tra using either the FICSB ell-derived or the FEBeCoP-derived 
beam agree very well (because the diff'erence between the two 
beams is small, as shown in Fig. 8). Figure 12 shows the results 



for all fields and frequencies. We also display in Fig. 13 the bin- 
bin correlation matrix, showing that two bins are not correlated 
by more than 10%. 

4.2.2. Systematic errors 

Our estimate of each power spectrum is aff'ected by diff'erent sys- 
tematic errors that must be accounted for separately from the sta- 
tistical errors derived in the previous section. These systematic 
uncertainties may introduce a bias in the final estimate and/or 
bias our Monte-Carlo estimate of the statistical uncertainties pre- 
sented above. We review here the diff'erent sources of these sys- 
tematics and evaluate their level. 



Mask impact- Our power spectrum estimation is performed on 
a limited sky patch. This induces power aliasing from angular 
scales larger than the size of the patch, an eff'ect that increases 
as the signal power spectrum steepens. POKER is designed to ac- 
count for this eff'ect (as well as the extra aliasing induced by 
holes in the map, if present) because the convolution kernel, 
Mhh', contains the information on mode coupling to the larger 
scales. We ran Poker on data maps that were embedded in larger 
regions that are zero-padded. There is no general prescription 
as to the size of the zero-padding that one should use, but for- 
tunately the results are very insensitive to the particular choice 
and whether or not the mask is apodized. By comparing different 
choices we found uncertainties at the 2% level, well below the 
statistical uncertainties. 



Template subtraction impact- Imperfect template subtraction 
will also lead to 'foreground' residuals that slightly bias our final 
estimate, P™, of the CIB anisotropy angular power spectrum. 
This residual level is given at first order by 



b' 



(29) 



with the pseudo- spectrum of the template and 6a the error 
on the global amplitude of this template. Figure 11 illustrates 
the level of these residuals for the particular case of the SP field. 
Although negligible compared to the statistical errors, they are 
accounted for in the error budget. 

Beam uncertainties- Uncertainty in the beam will also bias the 
estimate of the power spectrum because: 

1. The beam window function enters the computation of the 
convolution kernel M^^/ . Any beam error biases our estimate 
of Miyiy> and thus our final result. 

2. Beam uncertainties will translate into slight misestimation 
of A^™^(y) and C™^'^'(y), potentially biasing our final es- 
timate. 

Moreover, any beam uncertainty will also aff'ect our estima- 
tion of the statistical error bars. For example, the contribution 
of noise to the error bars scales roughly as (Ta^^ oc Nb/bl(y) (see 
noise curve on Fig. 11). As a consequence, any beam misesti- 
mation couples to the noise variance and leads to additional un- 
certainties on the final power spectrum estimate. This additional 
error is given at first order by 



Scr, 



Sh(y) 



bbiy) 



(30) 
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Figure 13. Modulus of the bin-to-bin correlation matrix derived from 
the simulation pipeline for the SP field at 353 GHz. Spatial frequency 
bins are not correlated by more than ^ 10%. 



where (Ta^^ is the noise error and Sbtiv) is the beam uncertainty. 
This may be important at small angular scales, depending on 
dbjjiy). A similar argument applies equally to sample variance. 

In this work we used the current, best determination of the 
beams. As explained in Sect. 3.2, the uncertainties on the scan- 
ning beam dominate these uncertainties. They are ~ 3% to ~ 6% 
of the FWHM, depending on the frequency (Planck HFI Core 
Team 2011a). 

We simulated the impact of such an error using the 
{simulation+analysis} pipeline presented in Sect. 4.2.1, assum- 
ing the appropriate frequency-dependent discrepancy between 
the beam used to simulate the maps and the beam used to anal- 
yse the maps. The bias induced by these uncertainties and their 
impact on the statistical error bars are derived by comparing the 
estimated power spectra and their MC-variance with and with- 
out the beam discrepancy. The bias so induced is the dominant 
uncertainty, larger even than the statistical error bars at small an- 
gular scales for measurements at 545 and 857 GHz. We took this 
bias into account in the modelling (Sect. 5.5). 

4.3. Robustness 

Many additional tests have been done to test the robustness of 
CIB anisotropy power spectra. First, instead of removing the 
CMB contribution and fitting only the Hi correlation coefiftcients, 
we searched for the best fit simultaneously using the Hi and 
CMB templates 

ivi^.y) = Yj KKM^y) +/^M^^y)cMB + Cy(x,y) . (3i) 

This allowed us to take into account the photometric inter- 
calibration uncertainties, which are about 2% for the CMB chan- 
nels (Planck HFI Core Team 2011a). We also fitted the low- 
frequency channels for only CMB 



Iy(x, y) = /32ly(x, y)cMB + Cy(x, y) . 



(32) 



Figure 12. CIB anisotropy power spectra of the six fields and their 
combined spectrum. 



We found that f5i and (52 agree at the -0.05% level. The 
diflTerence between the two coeflftcients and unity is less than 
1% (3%) at 217 GHz (353 GHz). The determination of the A is 
however very noisy because of the small area of our fields. They 
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are compatible with Pi^ = 1 at 217 GHz, though at 353 GHz 
they fall below unity by 2-3%. We did not correct for these 
inter-calibration coefficients and took conservative errors on 
the residual CMB contamination in our error pipeline (2% at 
217 GHz and 3% at 353 GHz, see Sect. 4.2). Fitting for both 
CMB and Hi or just Hi on CMB -cleaned maps changes the 
dust-Hi emissivities by less than 2%. This was again taken into 
account in the error simulation pipeline. 

To test the reliability of the noise power spectrum mea- 
surement, we recomputed the CIB anisotropy power spectrum 
using the cross -correlation between half-pointing period maps 
instead of the auto-correlation. We removed from each half map 
the same CMB and Hi (with emissivities taken as those of the 
total map). On average, for the range of i considered in this 
paper, the two methods agree at better than the 1% level (1% 
and 0.05% at 217 and 857 GHz, respectively). 

Finally, one way to asses the robustness of our determination 
is to compare the CIB anisotropy power spectra for our different 
fields. The fields all have different noise, cirrus, and CMB con- 
tributions and they are all independent. The comparison made in 
Fig. 12 shows that they are all compatible within error bars. 
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Figure 14. Field-combined CIB anisotropy power spectra at 217, 353, 
545, and 857 GHz. The dashed line shows the expected sum of the dusty 
and radio galaxy shot-noise power (from Table 3). The power spectra 
at 217, 353, and 545 GHz were arbitrary scaled to allow for a better 
comparison between frequencies (they were multiplied by 2x10^, 10^ 
and 10^, respectively). 



4.4. Power spectra from combined fields 

Our final estimate of the CIB anisotropy angular power spectra 
at different frequencies was derived by combining each power 
spectrum estimated on the six different fields 

C = Z<x^' (33) 

with A an index running over fields and an appropriate 
weight, to be defined below. The same binning was adopted for 
each field. The bin-bin covariance of the 'combined-field' es- 
timator, P^^'^^\ is a function of the bin-bin covariance of each 
'single-field' estimator, P^, and the covariance between two 
'single-field' estimators, and P^, follows 

C^^^ = Coy{pI''\p^?'^) (34) 

The error bars on P^^^^^ are simply given by -yjc^ and optimal 

weights are derived by searching for those minimizing these. 

Our MC simulations showed that different bins are not corre- 
lated by more than 10%, and our fields were widely enough sepa- 
rated that individual measurements are uncorrelated. Neglecting 
field- to-field and bin-to-bin covariances, the optimal weights be- 
come the usual inverse variance weighting: 

cr-l 

Wt = , (35) 

where crpA is the statistical error bars derived from the MC sim- 
ulations as previously described. 

The final CIB anisotropy power spectra estimates were there- 
fore computed by inserting the inverse variance weights in 



Eq. 33. The total statistical uncertainties were obtained by insert- 
ing those weights in Eq. 34, where Gov (p^, P^,) stands for the 
statistical uncertainties only (i.e. the standard quadratic summa- 
tion). The total systematic errors were obtained by linearly sum- 
ming the weighted systematic error on each field (any covariance 
between fields is neglected). We stress that irrespectively of the 
way weights are derived the full bin-bin covariance matrix could 
be computed^ . However, the forthcoming cosmological interpre- 
tation of the derived power spectra assumes a zero bin-to-bin 
correlation, and we only provide here the diagonal elements, i.e., 
error bars, of the final covariance. 

Our results are displayed in Fig. 14 and the data points are 
given in Table 4. Though our weighting is slightly suboptimal, 
the final angular power spectra are measured with high signal- 
to-noise compared to the single-field measurements. 

5. Overview and comparison with previous worl< 

The measured power and the shot noise predicted by the model 
discussed previously is shown in Fig. 14. Our measurements do 
not allow us to detect the shot-noise component, which will dom- 
inate on smaller scales than those probed here. However, the 
predictions are quite close to the measured power at the high- 
est indicating that further analysis of the CIB in Planck up to 
£ ~ 3000 might allow a measurement of the shot-noise compo- 
nent. The figure also reveals that the shape of the power spec- 
trum is remarkably similar at the four frequencies, being identi- 
cal within the 1 cr statistical errors for all data points but the last 
two at 217 GHz (i.e., £ ^ 1717 and 2060, which are 1.6 cr from 
the points at the other frequencies). This suggests that over the 
range of frequency and £ probed here the clustering properties 
do not evolve much and/or the galaxy populations responsible 
for CIB anisotropies are the same. We will return to this point 
in Sect. 5.5. We start in this section by analysing the frequency 
dependence of the CIB anisotropies and CIB mean, and compar- 

^ Assuming a vanishing bin-to-bin covariance for optimal weights 
computation does not prevent us from deriving the full bin-bin covari- 
ance matrix and just leads to sub-optimality in terms of error bars. 
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55.95 


43.65 


33.71 


23.18 


16.49 


16.67 


2.11 


20.51 


ACf^xlO^ 


28.23 


8.53 


6.08 


5.29 


4.90 


5.07 


5.59 


5.38 


6.39 


ACf xlO^ 


0.10 


0.18 


0.33 


0.46 


0.50 


0.52 


0.76 


0.13 


1.90 


353 GHz 




Q xlO^ 


755.17 


414.68 


215.18 


193.02 


117.35 


110.63 


92.76 


80.20 


74.84 


ACf ^ xlO^ 


83.18 


33.25 


16.05 


12.64 


8.64 


8.05 


6.84 


6.70 


7.56 


ACf xlO^ 


0.47 


1.04 


1.23 


2.01 


1.93 


2.65 


3.21 


3.99 


5.33 


545 GHz 




Q xlO' 


2246.54 


1091.71 


547.60 


454.01 


302.39 


271.87 


231.55 


196.63 


168.44 


ACf^xlO^ 


229.92 


81.11 


35.42 


24.90 


16.12 


13.15 


10.0 


8.12 


7.12 


ACf ^ xlO^ 


2.13 


4.16 


4.72 


7.14 


7.50 


9.84 


12.07 


14.76 


18.18 


857 GHz 




Ct xlO"^ 


561.65 


262.93 


120.57 


93.57 


67.65 


53.24 


46.82 


39.02 


33.88 


ACf ^ xlO"^ 


58.61 


18.82 


7.40 


4.79 


3.24 


2.36 


1.76 


1.35 


1.09 


ACf xlO"^ 


0.63 


1.18 


1.22 


1.72 


1.97 


2.26 


2.86 


3.44 


4.29 



Table 4. CIB anisotropy Q, at 217, 353, 545, and 857 GHz in jiK^y^^xsr. The conversion to Jy^sr ^ (with the photometric convention 
y/y=constant) involves multipHcation by 231483, 83135, 3391.5 and 4.99 at 217, 353, 545, and 857 GHz, respectively. ACf^ are the statisti- 
cal errors; AC^^^^ are the systematic errors introduced by the beam uncertainty (see Sect. 4.2.2). 



ing our anisotropy measurements with previous measurements 
at the same (or nearby) frequencies. 

5.1. Comparing the CIB mean and anisotropy SEDs 

The rms fluctuation in the CIB is related to the anisotropy power 
spectrum as 

(T^ = ^ f idiCe. (36) 
In J 

Table 5 gives approximations to this integral using the mea- 
sured values of over the range 200 < f < 2000. Statistical 
error bars on cr are computed with Monte Carlo simulations us- 
ing the statistical errors on the power spectra. The table also 
gives systematic errors (the second error term) corresponding 
to the photometric calibration uncertainties. Those values can 
be compared to the CIB absolute level. Cosmic infrared back- 
ground determinations based on FIRAS data from two groups 
can be used: 

- Fixsen et al. (1998) used three difl'erent methods to obtain 
the CIB mean. They average the three spectra to obtain one 
CIB mean spectrum, and then fit it by a modified black body. 
They find a dust temperature T = 18.5 ± 1.2 K, an opti- 
cal depth T = (1.3 ± 0.4) x 10"^ and an emissivity index 
J3 = 0.64 ± 0.12. The FIRAS spectrum is quite noisy so that 
the uncertainties on the parameters are large and the three 
parameters are highly degenerate. 

- Lagache et al. (1999) made two determinations of the CIB 
mean spectrum, using diff'erent methods to remove the fore- 
grounds (Hi and Hii) than Fixsen et al. (1998). One CIB 
mean spectrum is obtained in the Lockman Hole and one 
on 51% of the sky (to test isotropy). The two spectra agree 
very well (see Fig. 6 of Lagache et al. 1999). There is a re- 
finement of the measurement in Lagache et al. (2000), which 
agrees within errors with the previous measurement. Gispert 
et al. (2000) fit the Lockman Hole spectrum with a modi- 
fied black body to derive the CIB mean values and uncer- 
tainties at certain wavelengths (their Fig. 5). The best fit has 
T = 13.6 ± 1.5 K, T = (8.9 ± 2.9) x 10"^ andyS = 1.4 ± 0.2. 



These parameters are very diff'erent to those found by Fixsen 
et al. (1998), but because of their degeneracy, the spectrum 
is quite close to that of Fixsen et al. (1998) in the relevant 
frequency range (200-1500 GHz). 



We integrated the two CIB mean fits through the HFI 
bandpass filters to obtain the values for the absolute signal 
given in Table 5. The two determinations agree to better than 
10%, 1%, and 20% at 857, 545, and 353 GHz, respectively. At 
217 GHz they diff'er by 60%, but are compatible within errors. 
It is unknown which of the determinations is the best, because 
the errors on the spectrum are dominated by systematic eff'ects 
linked to foreground removal that are difficult to estimate. 
For the uncertainties listed in Table 5 we fixed T and ft to 
their best-fit values and considered only the uncertainty on r 
(since the errors on the three parameters are highly correlated). 
Comparison between the CIB mean and anisotropy SED is 
shown on Fig. 15. We see that the CIB anisotropy SED is well 
described by the CIB mean spectrum of Gispert et al. (2000) 
with the amplitude scaled by 0.15. The CIB mean spectrum of 
Fixsen et al. (1998) is flatter, with an 857/353 colour of 4.1 
compared to 5.3 and 5.4 and an 857/217 colour of 12, compared 
to 27 and 21, for CIB anisotropics and the Gispert et al. (2000) 
CIB mean, respectively. A steeper rise in the CIB anisotropy 
SED compared to the Fixsen et al. (1998) CIB mean has also 
been seen in Hall et al. (2010). These authors combine SPT and 
CIB anisotropy data from the literature and obtain an 857/217 
colour of ^ 25 at ^ = 3000. This compares very well to the 
CIB anisotropy colour obtained with Planck, integrated over 
200 < ^ < 2000. 



In conclusion, we do not see any evidence for a difl'erent CIB 
mean and anisotropy SED, which is consistent with the galaxies 
that dominate the CIB mean being those responsible for CIB 
anisotropics. The Gispert et al. (2000) fit of the Lagache et al. 
(1999) CIB mean spectrum (200 < v < 1500 GHz) describes 
both the CIB mean and the CIB anisotropy SED equally well. 
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217 GHz 


353 GHz 


545 GHz 


857 GHz 


C^obs 


(3.7 + 1.3 + 0.08) X 10--^ 


(1.9 + 0.3 + 0.04) X 10-' 


(5.9 + 0.8 + 0.40) X 10-' 


(1.0 + 0.1 +0.07) X 10-1 




(5.4 + 1.7) X 10"^ 


(1.6 + 0.5) X 10-1 


(3.7 + 1.1) X 10-1 


(6.5 + 2.0) X 10-1 




(3.4 + 1.1) X 10-' 


(1.3 + 0.4) X 10-1 


(3.7 + 1.2) X 10-1 


(7.1 +2.3) X 10-1 


CIBmod 


(3.2 + 0.6) X 10-' 


(1.2 + 0.1) X 10-1 


(3.5 + 0.2) X 10-1 


(6.3 + 0.3) X 10-1 



Table 5. RMS fluctuations in the CIB computed from Eq. 36 and CIB mean levels at 217, 353, 545 and 857 GHz. The subscripts 'obs' and 'mod' 
stand for observational and model values, respectively. The CIB model is taken from Bethermin et al. (2011). CIB^ and CIB§ are the Fixsen et al. 
(1998) and Gispert et al. (2000) best fits to the CIB spectra, respectively. The best fits and the model have been integrated in the HFI bandpasses. 
For the rms fluctuations both statistical and photometric calibration systematic errors are given. All numbers are given in MJy/sr for the photometric 
convention y/y= constant. 




Frequency (GHz) Multipole I 



Figure 15. Comparison of the observed CIB mean and anisotropy SED. 
The CIB measurements are from Lagache et al. (1999) (FIRAS spec- 
trum in black) and Penin et al. (2011b) (Spitzer and IRIS, pink dia- 
mond data points). The green and blue continuous (dashed) lines are the 
CIB fits from Gispert et al. (2000) and Fixsen et al. (1998) (multiplied 
by 0.15). The rms fluctuations of the CIB anisotropics, measured for 
200 < £ < 2000, are shown with the red dots. Their error bars include 
both statistical and photometric calibration systematic errors (linearly 
added), as given in Table 5. This figure shows that the CIB anisotropy 
SED is steeper than the Fixsen et al. (1998) best fit but very close to the 
Gispert et al. (2000) best fit. We see no evidence for diff'erent CIB mean 
and anisotropy SED. 



5.2. Comparison with SPT and ACT measurements 

SPT and ACT both measured the amplitude of CIB anisotropies, 
though at higher multipoles than those presented in this pa- 
per. The ACT measurement is on scales too small to be di- 
rectly compared with our measurements (the amplitude is given 
at ^ = 3000 while our last data point is at ^ = 2060). The 
SPT team computed the residual bandpowers after subtracting 
the tSZ, kSZ, CMB, and cirrus model components, quoting data 
points from £ = 2000 to 10,000. Fig. 16 shows the comparison of 
the HFI measurement at 217 GHz and the SPT measurements at 
220 GHz. Because the bandpass filters are not exactly the same, 
we applied a multiplicative correction factor (colour correction) 
to over-plot the SPT CIB anisotropy power spectrum on that of 
the HFI. This factor, the square of the HFI/SPT colour, is com- 
puted using the CIB SED of Gispert et al. (2000) convolved with 
the bandpass filters and is equal to 1.04. 

To interpret their data. Hall et al. (2010) used a phenomeno- 
logical model of CIB sources that assumed each galaxy has 
the same, non-evolving, modified blackbody SED and that their 



Figure 16. Comparison of SPT (Hall et al. 2010, dark open diamonds) 
and HFI measurements (red dots) at 217 GHz. The green dashed line 
corresponds to the SPT shot noise and the green dot-dashed line to the 
clustering model of Hall et al. (2010), the sum of the two is the continu- 
ous green line. The clustering model over-predicts by a factor ^ 2.4 the 
HFI power at ^ ~ 800. The blue dash-dotted line shows the clustering 
model divided by this factor. The clustering + shot noise (blue contin- 
uous line) now under-predicts the SPT data points, which may be the 
signature of non-linear contributions. 



light was a biased tracer of the mass fluctuations, calculated in 
linear perturbation theory. Moreover, the redshift distribution of 
the luminosity density was set by two parameters that fix the 
width and peak redshift. 

The green curve of Fig. 16 shows the Hall et al. (2010) 
model, normalized to the SPT bandpowers. We see that with 
this normalization, the power at large angular scales is larger 
by more than a factor of 2 than the HFI data. We also show as 
the blue curve the same model, except with amplitude adjusted 
to better agree with the HFI data. This downward adjustment of 
amplitude could arise from either a reduction in bias or in the 
amplitude of the mean CIB. This correction, of course, shifts the 
discrepancy to the smaller-scale SPT data. Since we expect the 
linear theory assumption will be better at large scales than at 
small, the discrepancy between model and data at small scales 
may be signaling the importance of non-linear corrections. 

5.3. Comparison with BLAST and SPIRE measurements 

Viero et al. (2009) presented BLAST power spectrum measure- 
ments at 1200, 857, and 600 GHz in the GOODS-South field. 
They detect CIB anisotropy and shot-noise power in the range 
940 < £ < 10,800. The measured correlations are well fitted by 
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a power-law over scales of 5-25', with A/// = 15.1% ± 1.7%. 
This level with respect to the CIB is the same as that found at 
the four HFI frequencies (see Sect. 5.1 and Fig. 15). Fitting to a 
linear theory power spectrum, they find that the BLAST galax- 
ies responsible for the CIB fluctuations have bias parameters, 
= 3.9 ± 0.6 and b = AA± 0.7 at 857 and 600 GHz, respectively. 
They further interpret their results using the halo model and find 
that the simplest prescription does not fit very well. One way 
to improve the fit is to increase the radius at which dark matter 
halos are truncated in the model (the virial radius) and thereby 
distribute satellite galaxies over a larger volume. They interpret 
this as being equivalent to having some star- forming galaxies at 
z > 1 located in the outskirts of groups and clusters. 

We show in Fig. 17 the comparison between the BLAST 
and HFI measurements at 857 and 545 GHz. Because the 
bandpass filters are quite diff'erent (particularly the 600 and 
545 GHz BLAST and HFI channels), we applied a colour 
correction as explained in Sect. 5.2, multiplying the BLAST 
CIB anisotropy power spectra by 0.7 and 1.05 at 545 and 
857 GHz, respectively. We see from Fig. 17 that the BLAST 
power spectra agree quite well with those from HFI, except that 
their largest- scale data points are systematically higher. This 
may be caused by contamination by residual Galactic cirrus 
emission in the BLAST power spectra. Also shown in the figure 
are shot-noise powers measured by BLAST. Once the colour 
corrections are applied, they are 1843 and 7326 Jy^ sr"^ at 545 
and 857 GHz, respectively. Their flux cuts are comparable to 
ours (they removed two sources above 400 mJy at 857 GHz, 
and no sources at 600 GHz). The measured shot noise levels 
are 1.6 and 1.2 times higher than the model predictions shown 
in Table 3 at 545 and 857 GHz, respectively. We also plot their 
best-fit halo model, which has a minimum halo mass required to 
host a galaxy of log(Minin/M0) = 11.5^^ ^ and an eff'ective bias 
Z^eff - 2.4. We see from Fig. 17 that their model is a very good 
fit to the HFI data points. Indeed, it provides a much better fit of 
the HFI data points than the BLAST data points ! 

We now compare our results with the HerschellSVl^^ mea- 
surements (Amblard et al. 2011). This comparison has to be 
made with caution, because the sources are masked down to a 
flux cut of 50 mJy in SPIRE, which is much smaller than the 
540 and 710 mJy flux cut used in HFI at 545 and 857 GHz, re- 
spectively. This large difl'erence in the source removal will afl'ect 
both the shot noise and the correlated components. We thus com- 
pare on Fig. 18 our HFI data points with a SPIRE measurement 
identical to the Amblard et al. (2011), but with no flux cut ap- 
plied (Amblard, private communication). In this figure, we only 
show the SPIRE measurements over the multipole range of 200 
to 3000 overlapping with Planck. With higher angular resolu- 
tion maps, SPIRE CIB anisotropy measurements extend down 
to sub-arcminute angular scales or ^ of ~ 2 x 10"^. For clarity 
we do not plot the small-scale power spectrum, but only concen- 
trate on the consistency between HFI and SPIRE CIB anisotropy 
at larger angular scales. We see from Fig. 18 that SPIRE mea- 
surements (green triangles) are significantly below the HFI CIB 
spectra (red dots). Two elements may explain this diff'erence: 

- Galactic cirrus: the cirrus signal in the SPIRE field is taken 
from existing measurements in the same field with IRAS 
100 //m and MIPS 160 yum, and the spectrum is extrapo- 
lated from 100 yum to SPIRE wavelengths using the spec- 
tral dependence of a Galactic dust model; this procedure 
is less accurate than the use of Hi, and overestimates the 
cirrus contamination because the IRAS data contain the 
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Figure 17. Comparison of BLAST and HFI measurements at 545 and 
857 GHz. HFI data points are the red circles; BLAST data points are 
the black diamonds. They were colour-corrected for the comparison 
(the colour was computed using the CIB SED of Gispert et al. (2000), 
integrated through the BLAST and HFI bandpass filters). The dashed 
line is the BLAST shot noise (also colour-corrected). Also shown is the 
BLAST best-fit clustering model (black dash-dotted line) and the total 
contribution (shot noise plus clustering; continuous green line). It pro- 
vides a good fit to the Planck data. Finally, we report in this figure a 
revised version of the SPIRE data points from Amblard et al. (2011) 
(blue triangles from Fig. 18, see text for more details). 



CIB anisotropies (Penin et al. 2011b). Indeed, we checked 
with our LH2 field, that overlaps with the SPIRE SWIRE- 
Lockman field to 43%, that in this region the cirrus contam- 
ination is negligible at the scales of interests for the com- 
parison (200 < £ < 2000). We accordingly added back the 
cirrus power spectra that were removed from the Amblard 
et al. (2011) power spectra. 
- HFI/SPIRE cross-calibration: SPIRE data are calibrated for 
point sources, with an accuracy of 15% (Swinyard et al. 
2010). The point-source to diflTuse-emission calibration con- 
version invokes an eflTective beam surface that is not perfectly 
determined yet. Planck/HFI is directly calibrated on diflfuse 
emission. As detailed in the Planck HFI Core Team (201 la), 
the accuracy is 7% at high frequencies. For now, Planck/HFI 
therefore has a more accurate photometric calibration for dif- 
fuse emission than SPIRE. At this stage, it appears that as- 
suming SPIRE beam surfaces corresponding to the integral 
of a Gaussian beam limited by diflfraction gives a more ac- 
curate SPIRE/HFI cross-calibration than the ofl&cial beam 
surfaces given in Swinyard et al. (2010). This preliminary 
cross-calibration has been established by comparing the dif- 
fuse emission from several Hershel surveys (some H ATLAS, 
SAG4, and Hi-Gal fields) to Planck/HFI data. Compared 
to the beam surfaces taken in Amblard et al. (2011)^, they 
will increase the power spectra by 10% and 20% at 857 and 
545 GHz, respectively. 

When corrected for the cirrus overestimate and the HFI/SPIRE 
cross-calibration on difiTuse emission, the SPIRE and HFI data 
points are now compatible (blue triangles and red dots on 



^ Amblard et al. (2011) use diff'erent values than those given in 
Swinyard et al. (2010), with beam surfaces of 1.77x10"^ sr and 
3.99x10"^ sr at 350 and 500 /um, respectively. 
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Figure 18. Comparison of SPIRE and HFI measurements at 545 and 
857 GHz in the overlapping multipole range. HFI data points are the red 
circles; SPIRE data points from Amblard et al. (201 1) are the black tri- 
angles. For these data points sources down to 50 mJy have been masked, 
there is consequently less power compared to HFI. The green triangles 
(Amblard, private communication) show the SPIRE CIB measurements 
identical to Amblard et al. (2011), but without a flux cut applied and 
thus they are directly comparable to the HFI measurement. They should 
agree with HFI, but are a factor ~1.7 and ~1.2 below the HFI CIB data 
points for 400 < ^ < 1500 at 857 and 545 GHZ, respectively. Indeed, 
they suffer from an overestimate of the cirrus contamination (by a factor 
2). Moreover, preliminary cross-calibration between SPIRE and HFI is 
increasing the Amblard et al. (2011) SPIRE power spectra by 10 and 
20% at 857 and 545 GHz, respectively (see Sect. 5.3 for more details). 
When corrected from these too factors (cirrus and cross-calibration), 
the SPIRE (blue triangles) and HFI measurements agree well. For this 
figure, all SPIRE data points were colour-corrected (colours were com- 
puted using the CIB SED of Gispert et al. (2000), integrated through 
the SPIRE and HFI bandpass filters). Error bars include only statistical 
errors (for SPIRE, error bars are only shown for the green triangles for 
sake of clarity). 



Fig. 18, respectively). The cirrus correction is the dominant ef- 
fect up to £ -500 and ^ -1500 at 545 and 857 GHz, respectively. 

5.4. A self-consistent, cosmological, IR, galaxy evolution 
model 

Our interpretation of the CIB anisotropy measurements with HFI 
relies on a model introduced in Penin et al. (2011a). The model 
builds upon the halo model formalism (see Cooray & Sheth 
2002, for a review) and populates dark matter halos with galax- 
ies using a HOD, modelling the emission of dusty galaxies us- 
ing the infrared evolution model of Bethermin et al. (2011, see 
Appendix C). Our main motivation for developing and using this 
parametric model is that it allows us to handle in a self-consistent 
manner the observational constraints coming from galaxy clus- 
tering and the CMB with more galaxy-evolution-centered mea- 
surements such as number counts or luminosity functions at vari- 
ous wavelengths and redshifts. This is a key feature of our model. 

Previous approaches, such as Amblard & Cooray (2007) and 
Viero et al. (2009), have used the Lagache et al. (2004) infrared- 
galaxy evolution model. Compared to Lagache et al. (2004) and 
Marsden et al. (2011), the parametric evolution of Bethermin 
et al. (2011) better reproduces the mid-IR to millimeter statisti- 
cal observations of infrared galaxies (number counts, luminosity 
functions, CIB, redshift distributions). This is important because 



we derive from this model the mean emissivity per comoving 
unit volume, introduced below, which is a key quantity for inter- 
preting CIB anisotropies. 

On the scales of interest to us we can use the Limber ap- 
proximation (Limber 1954) and write the angular (cross) power 
spectrum of infrared emission at two frequencies, v and y\ and 
at a multipole ^ as (e.g., Knox et al. 2001) 



dz 



X 



(37) 



where x is the comoving angular diameter distance to redshift 
z, a = (1 + z)~Ms the scale factor and jy(z) is the mean emis- 
sivity per comoving unit volume at frequency v and redshift z. 
The mean emissivity is derived using the empirical, parametric 
model of Bethermin et al. (2011)^: 



Jy(z) = (1 +Z) 



Jo 



dS S 



d^N 
dSdz 



(38) 



The remaining ingredient in the model is thus Pgg(k, z). As a 
foil to the HOD model for Pgg we begin with the simple, constant 
bias model in which 



Pgg(k,z) = bt^P,^(k,z), 



(39) 



where bn^ is a redshift- and scale-independent bias and Punik) 
is the linear theory, dark-matter power spectrum. We compute 
P\in(k) using the fit of Eisenstein & Hu (1998). We will see that 
this model is not sufl&cient to explain the CIB anisotropies that 
we measure. This is not unexpected; at the mean distance of the 
sources we are probing Mpc scales where non-linearities and 
scale-dependent bias are important. 

By contrast the HOD model computes P^^{k,z) as the sum 
of the contributions of galaxies within a single dark matter halo 
(Ih) and galaxies belonging to two diflTerent halos (2h): 



PUk) = Pi^(k) ^ P2^(k) . 



(40) 



The details of our assumptions for the Ih and 2h terms are given 
in Appendix C. On large scales P2h reduces to a constant bias 
(squared) times the linear theory power spectrum, while the 1- 
halo term becomes a scale-independent, shot-noise term. 

Before comparing our model to Planck observations, let us 
identify the parameters we hope to constrain with these data. 
The infrared galaxy evolution model of Bethermin et al. (2011) 
satisfyingly reproduces current number count observations and 
luminosity function measurements at the price of introducing 
a luminosity function characterised by thirteen parameters (see 
Appendix B). These thirteen parameters fully define the mean 
emissivities, jy(z), given in Eq. 38. The standard cosmological 
parameters (baryon density, tilt, etc.) mostly define the shape of 
the linear power spectrum in Eq. 39 and the geometric functions 
like x(z)- The HOD formalism we introduce in Appendix C re- 
quires four more parameters. Penin et al. (2011a) investigated 
this full parameter space and its degeneracies and concluded, 
not surprisingly, that the current generation of infrared galaxy 
clustering measurements will not allow us to constrain all these 
parameters simultaneously. Furthermore, they show that most of 
the constraints on the luminosity function evolution come from 
number counts and monochromatic luminosity function mea- 
surements. In the next section we threfore fix the luminosity 
function parameters to their best-fit values (from Bethermin et al. 
201 1) and vary only some of the HOD parameters. 



^ Note that for illustration purpose and where specified only, we will 
sometimes use the older phenomenological model of Lagache et al. 
(2004) (LDP). 
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5.5. Confronting the model with observations 

To confront the measurements with our model, we use a , — , 
Levenberg-Marquardt algorithm to perform a minimization. 
The;^^ for our data when compared to our model is given by 



Xy - ^ 

h 



model 



(/7)-Pf^(/7)> 



(41) 



Unless specified otherwise, we use errors including statistical 
and systematic photometric calibration errors (2, 2, 7 and 7% 
at 217, 353, 545 and 857 GHz, respectively, as defined in 
Sect. 2.1), added linearly, and we assume diagonal, uncorrected 
error bars as justified above. To reproduce the binning performed 
while measuring the power spectra, P^^^^^ih) is computed taking 
the average of £C^^^^^ at the minimum, maximum and mean 
i of each bin. We assume a Gaussian likelihood and assume 
that setting the fixed parameters (e.g., the luminosity function 
parameters) at their best-fit value is equivalent to marginalis- 
ing over them. It is important to note that the model power 
spectrum, P'^^'^^^ib) includes a shot-noise term (SN) defined in 
Sect. 3.1. Depending on the precise configuration we study, this 
shot-noise level is either fitted as an extra parameter or fixed to 
the predicted value. We also note that the models are derived 
for the Planck bandpass filters and are colour-corrected to be in 
Jy2 sr"^ (y/y =constant) as the data points, using the Gispert et al. 
(2000) CIB SED (colour corrections are (1.08)^ (1.08)^ (1.06)^ 
and (1.00)2 at 217, 353, 545 and 857 GHz, respectively). 

We remark here that this approach might be limited in two 
ways. First, we cannot exclude the possibility that our solution is 
a local extremum and not the global minimum of Eq. 41. While 
for all the results quoted above we checked the convergence with 
varying starting points, this particular point would have to be val- 
idated using for example a simulated annealing technique. We 
also did not explore the validity of the Gaussian likelihood ap- 
proximation. Some of these limitations could be resolved by im- 
plementing, for example, a Monte Carlo solver, and we defer this 
approach to future work. 

5.5.1 . Linear bias model and power-law constraints 

In this section we will illustrate the discussion looking at the 
545 GHz data only, but the same conclusions are reached with 
the other three frequencies. The relevant results are illustrated in 
Fig. 19. The data points correspond to our measurements with 
statistical error bars only. 

Fitting simultaneously for a shot-noise level and a constant 
bias, Z?iin, defined in Eq. 39, we obtain a good fit (x^/dof ^ 
0.36) as visible from the solid orange curve. The best fit bias is 
Z^iin = 2.45 ±0.18 which is consistent with previous results in the 
literature (Lagache et al. 2007; Viero et al. 2009). The cluster- 
ing contribution is plotted as the dashed orange line. However, 
the shot-noise level required to obtain this good fit (dot-dashed 
orange curve) is unrealistically high: (5.6 ± 0.7) x 10^ Jy^sr^. 
When compared to the expected level from our model whose 
68% C.L. amplitude is displayed as the yellow shaded area, our 
required level is ^ 5 times higher in power. Given that this model 
reproduces well all known number counts, the monochromatic 
luminosity function, CIB measurements (see e.g., the last line in 
Table 5), and Herschel/SFIRE shot-noise measurements (when 



^ The shot-noise levels required at the other frequencies to fit the 
data with the linear bias model are 27+6.7, 589+47, and 15915+1987 
Jy^sr"^ at 217, 353, and 857 GHz, respectively 
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Figure 19. In this plot we illustrate the constant bias model at 545 GHz. 
The orange (solid, dashed and dot-dashed) lines correspond to the best- 
fit linear model, its clustering component, and the shot-noise level, re- 
spectively. While this fit was performed using our fiducial emissivity 
defined in Eq. 38, for illustrative purpose we plot in green the analogous 
fit using the LDP emissivity. In both cases the required shot-noise level 
(dot-dashed lines) is well above the 68% C.L. predicted by Bethermin 
et al. (2011) and given in Table 5 (yellow contour). Conversely, the 
solid yellow line represents the best-fit curve when the shot-noise level 
is fixed to the expected value (Table 3) and only the constant bias is 
varied. The fit is obviously unsatisfactory. These results lead us to con- 
sider the linear bias model as unphysical, despite the good fit it provides 
(X^/dof ^ 0.36 (2.52/7)). For illustration purpose, we also show our 
best-fit power-law model defined in Eq. 42 (blue solid line). 



Frequency 


A 


n 


Reduced 


(GHz) 


Jy^ sr~^ 




(x'/dof) 


217 


51 + 5 


-1.04 + 0.13 


0.98 (6.92/7) 


353 


1117 + 46 


-1.03 + 0.06 


0.86 (6.07/7) 


545 


(114 + 7)xl02 


-1.09 + 0.10 


0.21 (1.51/7) 


857 


(35 + 2)xl03 


-1.18 + 0.10 


0.25 (1.75/7) 



Table 6. Power-law model best-fit parameters for each frequency as 
well as the reduced The errors corresponds to the Icr Gaussian er- 
ror including statistical and photometric calibration systematic contri- 
butions. 



no sources are removed, see Sect. 5.3), this excess level is ex- 
cluded (see also the discussion in Sect. 3.1). We thus conclude 
that the linear bias model is not a physically realistic fit to our 
data. The reason for this will be made clear in the next section, as 
it will appear clearly that the shot-noise level we fit here absorbs 
a strong non-linear component. 

To further illustrate this point, we show the solid yellow 
curve that represents the best-fit model when we fix the shot 
noise to the level predicted by our model and only vary bn^. 
The fit is obviously unsatisfactory. We also note that the spe- 
cific emissivity density that comes from the underlying model 
is unlikely to be the culprit. The solid, dashed, and dot-dashed 
green curves are the analogous fit when we use the emissivity 
coming from LDP. The conclusions remain unchanged. 

Galaxy correlation functions can be reasonably well fitted, 
over a limited range of scales, by power-laws. A power-law cor- 
relation function would project into a power-law angular corre- 
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lation function, so we also consider a power-law fit to our data, 

This simple two-parameter model (A,n) is a reasonable fit at 
all frequencies, giving a reduced ;^^/ Jo/ ^ 0.21 (1.51/7) at 
545 GHz. The best-fit values are given in Table 6 and the best-fit 
model at 545 GHz is displayed as the blue solid line in Fig. 19. 
Although the power-law model provides a good fit to the data, it 
does not provide physical insight into the properties of the galax- 
ies it describes. 

5.5.2. HOD model constraints 

We now consider the HOD model introduced in Sect. 5.4. In the 
most general configuration, our parametrisation of this model 
allows for four diff'erent parameters: Mmm, Msat, o^sat and (XiogM 
(see Appendix C). The full exploration of this parameter space 
turns out to be a difficult task beyond our scope in this paper. 
Therefore, we will restrict ourselves to only two parameters, 
Mmin and c^sat, which describe the mass above which we expect a 
halo to host a CIB -contributing galaxy and the slope of the high- 
mass end of the HOD. By varying these parameters we control 
the mean, galaxy-weighted halo mass and the satellite fraction 
in the model, which in turn control the amplitudes of the Ih and 
2h terms. 

We impose Msat = 3.3 Mmin (Penin et al. 2011a) and choose 
o'XogM = 0.65, motivated by the clustering observed in optical 
surveys (see Tinker & Wetzel 2010, for discussion and refer- 
ences). We did not find (XiogM to be a critical parameters for our 
fit, but letting it vary drives the fit to physically unrealistic region 
of parameter space. For each frequency, we fixed the Poisson 
noise level to the one expected from our model (see Sect. 3.1). 
To add to the robustness of our interpretation, we took one more 
conservative step. 

As stated above, our interpretation of these data requires 
the knowledge of the emissivity. While our fiducial model is 
well tested and reproduces all relevant current observations 
(Bethermin et al. 2011), these do not extend beyond a redshift 
of about 3.5. As such, extrapolations to higher redshifts are un- 
constrained by previous observations except for the integral con- 
straints provided by the CIB measurement discussed in Sect. 5.1. 
To let our conclusions be as model-independent as possible and 
also to isolate and constrain the high-z contribution, we made the 
extra assumption that the emissivity, j, is constant for z > 3.5 and 
we fitted for it simultaneously while solving for logjo Mmm and 
a sat- More precisely, we rewrite Eq. 37 as 

f 3.5 ^2 _ 

Jo dz X 

+ (jeffT £ *f P^H^^^ = ^> ' ^43) 

where we have introduced the eff'ective redshift-independent 
emissivity, j^^, that we will also solve for. We adopt the arbi- 
trary but reasonable z = 7 cut-off" of Bethermin et al. (201 1). 

We treated the four frequencies as independent and per- 
formed a single minimization per frequency. The results are il- 
lustrated in Fig. 20. The solid orange line represents the best- 
fit model per frequency. Our three-parameter model obviously 
fits each frequency very well. The orange dashed line represents 
the 2-halo (linear) term, while the orange dot-dashed line repre- 
sents the 1-halo (non-linear) term. The green curve corresponds 



to the assumed Poisson noise level. Clearly, the angular scales 
we probe require a modelling of both the linear and the non- 
linear contribution to the power spectrum for all frequencies, and 
the 1-halo term is similar in slope to the shot-noise term, which 
leads to a model degeneracy. 

Quantitative results are given in Table 7, where we quote a 
reduced as a goodness-of-fit measure. The errors quoted in 
Table 7 correspond to the Gaussian errors computed from the 
Fisher matrix at the best-fit values. Each I bin at any given fre- 
quency is considered independent from the others at all frequen- 
cies. For reference, we also give the results of a fit where we 
fixed i to the value given by the Bethermin et al. (2011) model 
and fitted for only logjg Mmin and Ofsat- While the best-fit values 
are consistent between the two models, it is clear from the ta- 
ble that allowing jeff to vary degrades strongly the constraints on 
log 10 Mmin and Ofsat- In fact, a strong degeneracy between y'eff and 
logio ^min is observed, as might be expected. 

As is the case with optical data, our data do not appear to 
require a departure from asat=l- We considered diff'erent ratios 
of Msat/Mmin, 2, 5, 10 or 20. None of them provided a better 
fit, and most of them required similar values of cKsat (see Tinker 
& Wetzel 2010, for a recent summary). 

Taking our emissivity model at face value, the best-fit angu- 
lar power spectrum at ^ = 2000 receives the following redshift 
contribution at 217 (353/545/857) GHz: 5% (4/34/71) between 
< z < 1, 7% (7/23/22) between 1 < z < 2, and 88% (89/43/7) 
for 2 < z. 

It is obvious from our Fig. 20 that the non-linear contribu- 
tion is degenerate with the Poisson noise level. This explains the 
problem faced by the linear model discussed in Sect. 5.5.1. Our 
data by themselves are not sufficient to explore this degeneracy 
and we thus rely on our model. For similar reasons, we do not 
discuss details of the implementation of the Ih term (e.g., the 
truncation radius in m(^, M)) unlike Viero et al. (2009). We ex- 
pect that a future joint analysis with higher angular resolution 
measurements from Herschel and SPT/ACT will allow us to al- 
leviate this degeneracy. 

As described in Sect. 5.1, our fiducial model predicts a mean 
CIB consistent with observations. While it is not possible to 
translate our constraints on (a weighted integral of f over 
redshift) into a prediction for the integral of j with the diff'erent 
weighting required to compute the CIB mean, rough estimates 
suggest that our values are consistent with the FIRAS measure- 
ments. 

The relatively good consistency between the best-fit values 
of log 10 ^min and CKsat obscrvcd across frequencies raises an in- 
teresting question: does a single model fit all our data? Or to 
put it another way, are the diff'erences between each frequency 
HOD subsantial? Diff'erent frequencies, loosely speaking, corre- 
spond to diff'erent redshifts for the dominant galaxy population. 
As such, consistency in log^g ^min and Ofsat could imply that the 
CIB fluctuations arise from a single subset of galaxies whose 
redshift evolution we capture well with our emissivity model, 
mass function, and HOD prescription, which is constant in red- 
shift. To illustrate this hypothesis. Fig. 21 shows for each fre- 
quency's best-fit model its prediction at 545 GHz. Even though 
this plot does not convey the uncertainty associated with each 
prediction, it clearly illustrates that each HOD leads to substan- 
tially diff'erent predictions and thus that the frequency depen- 
dence of the HOD may be significant. We postpone to future 
work more quantitative statements on the implications for the 
clustering of galaxies at high redshift. 
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Figure 20. Angular power spectrum of CIB anisotropies at 217, 353, 545, and 857 GHz. Each panel corresponds to one frequency. For each 
frequency, the blue points correspond to the angular auto power- spectra, and the associated error bars include statistical and photometric calibration 
systematic contributions. The best-fit model per frequency (including shot noise) corresponds to the solid orange line. The dashed (dot-dashed) 
orange lines correspond to the 2h (Ih) contributions. The green triple dot-dashed curve corresponds to the Poisson noise level, fixed to its expected 
value. To obtain these fits, three parameters per frequency were varied: log^o ^min' ^sat and jeff • The fits are obviously qualitatively very good. 



Frequency (GHz) 


logioMmin [/?"^Mo] 




Jeff [Jy/Mpc/sr] 


Reduced (X^/dof) 


217 


11.95 + 2.10 


1.30+1.16 


7.51 +0.75 X 10^ 


2.68(16.1/6) 


353 


12.49 + 0.42 


1.39 + 0.42 


2.00 +0.29 X 102 


2.42 (14.5/6) 


545 


12.35 + 1.01 


1.17 + 0.65 


3.11 +3.85 X 10^ 


0.50 (3.04/6) 


857 


12.20 + 0.51 


1.02 + 0.87 


3.14 +17.0 X 10^ 


0.73 (4.40/6) 


217 


11.82+ 1.92 


1.17 + 2.38 


N/A 


1.14(7.96/7) 


353 


12.50 + 0.09 


1.35 + 0.20 


N/A 


0.80 (5.64/7) 


545 


12.35 + 0.94 


1.17 + 0.45 


N/A 


0.35 (2.46/7) 


857 


12.21 + 1.23 


0.96 + 0.73 


N/A 


0.60 (4.22/7) 



Table 7. Best-fit values for each frequency, as well as the reduced ;^2. The errors correspond to the Icr Gaussian errors, including statistical and 
photometric calibration systematic contributions. Systematic errors introduced by the beam uncertainty (see Sect. 4.2.2) are not included here, but 
contribute less than an extra 10% to the error budget. The upper half of the array allows for a freely varying j^^ per frequency, while in the bottom 
half Jeff is fixed to the extrapolation coming from our model. 



6. Conclusion 



We presented the first measurement of CIB anisotropies with 
Planck, detecting power from 10' to 2°. Owing to the excep- 
tional quality of the data, and using a complete analysis of the 
diflTerent steps that lead to the CIB anisotropy power spectra, we 



were able to measure the clustering of dusty, star- forming galax- 
ies at 217, 353, 545, and 857 GHz with unprecedented precision. 

We worked on six independent fields, chosen to have high 
angular-resolution Hi data and low foreground contamination. 
The CIB maps were cleaned using templates: Hi for Galactic 
cirrus; and the Planck 143 GHz maps for CMB. Having Hi data 
is necessary to cleanly separate CIB and cirrus fluctuations. 
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Figure 21. Predicted 545 GHz power spectra derived from each fre- 
quency's best-fit model. For the best-fit model at 217, 353, 545 and 
857 GHz, we plot the predicted clustering plus shot-noise power spec- 
tra at 545 GHz. Also shown are the HFI data points at 545 GHz (red 
diamonds). This plot suggests that the fits across frequencies are fairly 
diff'erent, which hints at an evolution in the population of galaxies we 
probed. We note, however, that the uncertainties associated with each 
prediction are not fully characterised by our method. 



Because the CIB anisotropies and Galactic dust have similar 
SEDs, blind component separation methods do not adequately 
distinguish CIB anisotropies from cirrus emission. The 143 GHz 
Planck channel, cleaned from sources and filtered, provides a 
good template for the CMB because it has instrument noise 
and an angular resolution close to the higher frequency channels 
from which we measure the CIB. It also has the advantage of 
being an"internal" template, meaning its noise, data reduction 
process, photometric calibration, and beam are well known. 

We obtained CIB anisotropy maps that reveal structures pro- 
duced by the cumulative emission of high-redshift, dusty, star- 
forming galaxies. The maps are highly correlated at high fre- 
quencies. They decorrelate at lower frequency, as expected from 
models of the redshift distribution of sources producing the 
CIB anisotropies (e.g. Fernandez-Conde et al. 2008; Penin et al. 
2011a). In these models, at 217 GHz the contribution of z > 2 
galaxies is becoming dominant, while at higher frequencies the 
dominant sources are at lower z. We computed the power spec- 
tra of the maps and their associated errors using a dedicated 
pipeline, based on the POKER algorithm (Ponthieu et al. 2011). 
After a careful examination of many systematic effects, use of 
the best determination of the beam window function and instru- 
ment noise determined from jack-knife methods, we ended up 
with measurements of the angular power spectrum of the CIB 
anisotropy, Q, at 217, 353, 545, and 857 GHz, with high signal- 
to-noise ratio over the range 200 < i < 2000 (see Table 4 and 
Fig. 12). 

The SED of CIB anisotropies is not diflTerent from the CIB 
mean SED, even at 217 GHz. This is expected from the model of 
Bethermin et al. (2011) and reflects the fact that the CIB mean 
and anisotropies are produced by the same population of sources. 
Our measurement compares very well with previous measure- 
ments at higher at 220 GHz by SPT (Hall et al. 2010) and at 
600 and 857 GHz by BLAST (Viero et al. 2009). On the contrary, 
HerschellSVl^^ measurements from Amblard et al. (2011), but 
with no sources removal to be comparable to HFI (Amblard, pri- 
vate communication), are significantly lower than our measure- 



ments at high frequencies owing to an overestimate of the cirrus 
contribution. Preliminary cross-calibration between SPIRE and 
HFI is also increasing the Amblard et al. (2011) power spectra 
by 10 and 20% at 857 and 545 GHz, respectively. 

From the Planck data alone we can exclude a model where 
galaxies trace the linear theory matter power spectrum with a 
scale-independent bias: that model requires an unrealistic high 
level of shot noise to match the small-scale power we observe. 
Consequently, we developed an alternative model that couples 
the dusty galaxy, parametric evolution model of Bethermin et al. 
(2011) with a halo model approach. Characterised by only two 
parameters, this model provides an excellent fit to our measured 
anisotropy angular power spectrum for each frequency treated 
independently. Whereas in principle these two parameters could 
oflTer us unique insights into the clustering and the nature of 
dusty galaxies at high redshift, the current uncertainties in the 
underlying model prevent us from drawing detailed inferences. 
Our results suggest that a diflTerent HOD is required at each fre- 
quency, which is consistent with the fact that we expect each 
frequency to be dominated by contributions from diflTerent red- 
shifts. We find that half of the contribution to the power spec- 
trum at ^=2000 comes from redshifts lower than 0.8 and 1.5 at 
857 and 545 GHz, respectively. Those numbers are quite robust 
against exact evolution of dusty galaxies comoving emissivity at 
high-redshift {z >3.5). This is not the case at lower frequencies 
and our best-fit model predicts that about 90% of the anisotropies 
power at ^=2000 come from redshifts z >2 at 353 GHz and 
217 GHz. 

Further modelling and interpretation of the CIB anisotropy 
will be aided by the use of cross-power spectra between bands, 
and by the combination of the Planck and Herschel data at 857 
and 545/600 GHz and Planck and SPT/ACT data at 220 GHz. 
This combination will measure the CIB anisotropy power spec- 
trum over a wide range of scales, covering the three regimes 
where we expect the 2-halo, 1-halo and shot-noise contributions 
to dominate. More progress could be made by measuring the 
CIB anisotropies over more sky and at lower frequencies (at 
least 143 GHz) with Planck. Going to lower frequency extends 
our reach in redshift, and is also important for CMB analysis 
and measurement of the SZ power spectrum. Additional infor- 
mation will be obtained by cross-correlating the CIB maps with 
external tracers of the density field, like the galaxy and quasar 
distributions in large area catalogues (such as those from the 
SDSS, VIKINGA^ISTA and KIDS/VST surveys). This will ad- 
ditionally constrain 1) the populations contributing most to the 
CIB; and 2) the relative bias between the external tracer and 
the distribution of far-infrared emission. A particularly interest- 
ing cross-correlation may be between the CMB lensing conver- 
gence and the CIB maps (Song et al. 2003). The lensing and CIB 
anisotropies are expected to have a high degree of overlap, and 
lead to a signal readily detectable by Planck. This signal will 
give a direct and independent measure of the bias. 
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Appendix A: From HFI maps to CIB power spectra: 
flow charts of the different steps 

We summarize with the two flow charts presented in Fig. B.l and 
B.2 the procedures of data preparation and cleaning, and power 
spectra measurements and errors evaluation. 

Appendix B: The parametric backward evolution 
model of dusty star-forming galaxies 

In this appendix we give some details about the dusty star- 
forming galaxies evolution model we are using in our modelling 
of CIB anisotropics. The model is fully described in Bethermin 
et al. (2011). Two ingredients come into play: the spectral en- 
ergy distribution (SED) of galaxies and the luminosity function 
(LF) evolution. 

The SEDs are from the templates library of Lagache et al. 
(2004), which consists of two galaxy populations: a star- 
forming galaxy population with SEDs that vary with IR 
bolometric luminosities, and a normal-galaxy population with a 
template SED that is fixed and is colder than the star-forming 
galaxy templates and is scaled with IR bolometric luminosities. 
The normal and star-forming galaxies are dominant at low- 
and high-luminosity, respectively. The fraction of each galaxy 
population as a function of the bolometric luminosity was given 
using a smooth function 



starburst 
~0 



1 -h ^/l[logio(LiR/Lpop)/crpop] 



(B.l) 



where th is the hyperbolic tangent function, Lpop the luminosity 
at which the number of normal and star-forming galaxies are 
equal, and (Xpop characterises the width of the transition between 
the two populations. At Ljr 
50%. 



-pop 5 



the starbursts fraction is 



They assume that the luminosity (LF) is a classical double 
exponential function: 



O(Lir) = 



x(^)^- 



(B.2) 



where O(Lir) is the number of sources per logarithm of lumi- 
nosity and per comoving volume unit for an infrared bolometric 
luminosity Ljr, is the normalization constant characterising 
the density of sources, is the characteristic luminosity at the 
break, and \ - a and 1 - a - I / cr^ /Iri^ (10) are the slope of the 
asymptotic power-law behaviour at low and high luminosity 
respectively. 
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Figure B.l. Cartoon illustrating the different steps from HFI frequency 
maps to CIB power spectra for one field. 



A continuous LF redshift-evolution in luminosity and density is 
assumed following oc (1 + zY^ and oc (1 + zY'^, where ri 
and r^p are parameters driving the evolution in luminosity and 
density, respectively. It is impossible to reproduce the evolution 
of the LF with constant and r^. They consequently authorize 
their value to change at two specific redshifts. The position 
of the first redshift break is a free parameter and converges to 
the same final value (z ~ 0.9) for initial values < z < 2. To 
avoid a divergence at high redshift, the second break is fixed 
at z=2. The position of the breaks are the same for both tl and r^. 

The model has 13 free parameters that are summarized 
in Table C.l. The parameters were determined by fitting the 
model to published measurements of galaxy number counts and 
monochromatic LF measured at given redshifts: 

- Number counts: Spitzer counts at 24, 70 and 160 yum, 
Herschel counts at 250, 350 and 500 yum, and AzTEC counts 
at 1.1 mm. 

- Monochromatic LFs: IRAS local LF at 60 yum, Spitzer LF at 
24 yum at z=0, at 15 yum at z=0.6, at 12 yum at z=l, and at 
8 yum at z=2. 

- FIR AS CIB spectrum between 200 jim and 2 mm. 

Measured redshift distributions were not used because the 
cosmic variance and the selection effects were currently poorly 
quantified. The best-fit parameters as well as their uncertain- 
ties and degeneracies were obtained using a Monte Carlo 
Markov chain (MCMC) Metropolis-Hastings algorithm. The 
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Figure B.2. Cartoon illustrating the angular power spectrum measurement and errors estimate using the POKER algorithm (Ponthieu et al. 2011). 



model adjusted on deep counts and monochromatic LFs at 
key wavelengths also reproduces recent very discriminating 
observations well, such as the Jauzac et al. (2011) measured 
redshift distribution of the CIB. 

We used the so-called mean model, which is ob- 
tained using the mean value of the parameters as 
given in Table C.l without the lensing contribution 
(dndsnudz_arr_nolensing_meanmodel -final file on the 
http://www.ias.u-psud.fr/irgalaxies/ web page). 



Appendix C: The halo model 

In this appendix we give the details of our halo modelling. 
Neglecting scale-dependent halo bias, the distinction between 
central and satellite galaxies and halo exclusion, and assuming a 
Poisson distribution of galaxies, the Ih and 2h terms in Pgg have 
simple analytic expressions (Cooray & Sheth 2002): 



dN (A^gal(A^gal - 1)) 

dM , 



u^(k, M); 



[J.M^KM) 



(A^gal) 



u(k, M) 



(C.l) 
(C.2) 



Here M is the halo mass, dN/dM is the halo mass function, 
u(k, M) is the Fourier transform of the (normalized) halo den- 
sity profile, b{M) the halo bias and (A^gai) is the mean number of 
galaxies in a halo of mass M. 

The mean number density of galaxies, ^gai, can be written 



^gal 



^M^<A.g,). 



(C.3) 



Note that on large scales uik ^ 0, M) ^ 1 so that we can define 
the "eff'ective" bias as 



f dM —b(M) 
J dM ^ ' 



(A^gal) 
^gal 



(C.4) 



and the 1-halo term becomes scale-independent {i.e., a shot- 
noise term). 

We used the fitting function of Eisenstein & Hu (1998) to 
compute Piin, an NFW profile (Navarro et al. 1997) truncated at 
the virial radius to compute uik, M), and we relied on the mass 
function fit of Tinker et al. (2008) with its associated halo bias 
prescription (Tinker et al. 2010). All these relations were cali- 
brated through the use of N-body simulations. Our definition of 
halo mass is the mass interior to a radius within which the mean 
density is 200 times the mean density of the Universe. 
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Parameter 



Description 



Value 



a 
cr 

L^(z=0) (xlQio Lo) 

(f>^ (z=0) (xlO-3 gal/dex/Mpc^) 

^phii,,lz 
Zbreak,l 

^phii, ,mz 
Zbreak,! 

(XIO'O Le) 



Faint-end slope of the infrared bolometric LF 1.223 + 0.044 

Parameter driving the bright-end slope of the LF 0.406 + 0.019 

Local characteristic luminosity of the LF 2.377 + 0.363 

Local characteristic density of the LF 3.234 + 0.266 

Evolution of the characteristic luminosity between and Zbreak,i 2.93 1 + 0. 1 19 

Evolution of the characteristic density between and Zbreak,\ 0.174 ± 0.196 

Redshift of the first break ' 0.879 + 0.052 

Evolution of the characteristic luminosity between Zbreak,i and Zbreak,! 4.737 ± 0.301 

Evolution of the characteristic density of between Zbreak,i and Zbreaka -6.246 + 0.458 

Redshift of the second break ' ' 2.000 (fixed) 

Evolution of the characteristic luminosity for z> Zbreaka 0.145 + 0.460 

Evolution of the characteristic density for z> Zbreak,i -0.919 + 0.651 

Luminosity of the transition between normal and starburst templates 23.677 + 2.704 

Width of the transition between normal and starburst templates 0.572 + 0.056 



Table C.l. Dusty star-forming galaxies evolution model parameters, fitted to selected infrared observations (table extracted from Bethermin et al. 
20 11). The errors are derived from a MCMC analysis. 



The HOD describes the way galaxies populate the dark mat- 
ter halos. While we do not distinguish between central and satel- 
lite galaxies in the above, the functional form we adopt for the 
mean occupation is modelled on the form frequently used in op- 
tical observations (e.g., Zheng et al. 2005) 



with 



1 -Kerf 



log M - log Mn: 



and 



1 -h erf 



logM-log2M^ 



(C.5) 



(C.6) 



(C.7) 



These definitions ensure that M^in is the halo mass at which a 
halo has a probability of 50% of having a central galaxy. We in- 
troduce (TiogM to allow scatter in this relation between halo mass 
and observable, which is important on large scales. Following 
Zheng et al. (2005), we assume a Poisson distribution for A^^^^ 
and write 



<A^gal(A^gal-l)> = 2<7Vsat> + <A^sat>' 



(C.8) 



Within this parametrisation, halos with M «c Mmin will not host 
any galaxies, whereas those with M » Mmin are almost certain 
to contain one. The satellite occupation has a similar cut-oflf, but 
the mass is chosen to be twice Mmin, so that halos with a low 
probability of having a central galaxy are unlikely to contain a 
satellite galaxy (see Tinker & Wetzel 2010, for a further discus- 
sion of this form). 

We note that this parametrisation was introduced to repro- 
duce the observed clustering of luminosity-threshold samples of 
optical galaxies at < z < 2. We are therefore making substan- 
tial assumptions when applying this same parametric form to 
dusty star-forming galaxies at higher redshift. As we shall see, 
however, our constraints on even this form of the HOD are weak 
enough to argue against introducing additional degrees of free- 
dom in the model. 
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